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ABSTRACT 

We study the role of ambipolar diffusion (AD) on the non-linear evolution of the MRI in protoplan- 
etary disks using the strong coupling limit, which applies when the electron recombination time is 
much shorter than the orbital time. The effect of AD in this limit is characterized by the dimension- 
less number Am, the frequency of which neutral particles collide with ions normalized to the orbital 
frequency. We perform three-dimensional unstratified shearing-box simulations of the MRI over a 
wide range of Am as well as different magnetic field strengths and geometries. The saturation level 
of the MRI turbulence depends on the magnetic geometry and increases with the net magnetic flux. 
There is an upper limit to the net flux for sustained turbulence, corresponding to the requirement 
that the most unstable vertical wavelength be less than the disk scale height. Correspondingly, at 
a given Am, there exists a maximum value of the turbulent stress amax- For Am < 1, the largest 
stress is associated with a field geometry that has both net vertical and toroidal flux. In this case, 
we confirm the results of linear analyses that show the fastest growing mode has a non-zero radial 
wave number with growth rate exceeding the pure vertical field case. We find there is a very tight 



l/2a at the saturated 



correlation between the turbulent stress a and the plasma 

state of the MRI turbulence regardless of field geometry, and amax rapidly decreases with decreasing 
Am. In particular, we find amax ~ 7 x 10"'^ for Am = 1 and amax ~ 6 x 10^^ for Am = 0.1. 
Subject headings: magnctohydrodynamics — instabilities — methods: numerical — planetary systems: 
protoplanctary disks — turbulence 



1. INTRODUCTION 

One of the most fundamental questions in accretion 
disk dynamics is how the disk transports angular mo- 
mentum and accretes to the c entral object. The magn e- 
torotational instability (MRI, iBalbus fc Hawlevlll991| ) is 
widely considered to be the most likely mechanism for the 
transport process. The non-linear evolution of the MRI 
under ideal MHD conditions has been studied extensively 
using both local (iHawlev et al.iri 99 5: St one et al. 19961: 
Miller & Stone 200 0) and globS7Armitagd ll998l:lHawlevl 
200a i2001c .Fromang fc NelsonI 120061 ) numerical simula- 
tions. It has been found that MRI generates vigorous 
MHD turbulence and produce efficient outward transport 
of angular momentum whose rate is compatible with ob- 
servations. However, accretion disks in some systems are 
only partially ionized, and non-ideal MHD effects need 
to be taken into account. In particular, most regions 
of the protoplanctary disks (PPDs) are too cold for suf- 
ficient thermal ionization, and effective ionization may 
be achieved only in the disk surface layer due to exter- 
nal ionization sources such as X-ray rad iation from the 
central star and cosmic ray ionization (jCammi ^ [19961 : 
iStone et all l2000f ). Non-ideal MHD effects reflect the 
incomplete coupling between the disk material and the 
magnetic field, and substantially affect the growth and 
saturation of the MRI. 

There are three non-ideal MHD effects as manifested in 
the generalized Ohms's law, namely the Ohmic resistiv- 
ity. Hall effect and ambipolar diffusion (AD), with three 
different regimes associated with the relative importance 
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of these terms. In general, the Ohmic term dominates 
at high density and very low ionization, the AD term 
dominates in the opposite limit, and the Hall term is 
important in between. So far the majority of studies 
have been focused on the Ohmic regime. In this case, 
MRI is damped when th e Elsasser number A = y-j^/rir t 
falls below order unity ([jh] 119961 : iTurner et ai1 |2007^. 
where vaz is the Alfven velocity in the vertical direction, 
77 is the Ohmic resistivity, il is the angular frequency 
of Keplerian rotation. Another often quoted criterion is 
the magnetic Reynolds number Re a/ = cl/rjil, > 10"* for 
MRI to be self-sustained (where Cg is the sound speed), 
which has the advantage of being independent of the 
magnetic field (jFleming et al.|[2000| ). Ohmic resistivity 
has been used extensively to model the layered accre- 
tion in PPDs, where the surface layer of the disk is 
sufficiently ionized to couple to the magnetic field and 
drive the MRI turbulence, while the midplane region is 
too poorly ionized and "dead" (iFleming fc Stond 120031: 



Turner et al .1120071: [Til^r fc Sanol l2008l: rilgner fc NelsonI 
200a lOishi fc Mac Lowll2009l) ! 

The importance of Hall and AD terms in PPDs has 
been studied in a number of theoretical works, but rel- 
atively little attention has been paid to numerical simu- 
lations of the non-linear regime. Linear analysi s of the 
MRI i n the Hall regime have b e en pe rformed by IWardld 
([1999}) and IBalbus fc Tergueml (|200l . The growth rate 
is strongly affected by the Hall term and depends on 
the sign of B • il. Nevertheless, numerical simulations 
including both the Hall and Ohmic terms (where the 
Ohmic term dominates) showed that the Hall term does 
not strongly affect the saturation amplitude of MRI 
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(jSano fc Stonell2002allb[ ). It is yet to study the behavior 
of MRI in the regime where HaU effect dominates over 
other terms, and to include the HaU term in the more 
rcahstic vcrticaUy stratified simulations. 

The relative motion between ions and neutrals leads to 
AD. AD is ideally studied using the two-fluid approach, 
where the ions and neutrals are treated as separate flu- 
ids, coupled by the ion- neutral drag via collisions. More- 
over, ion and neutral densities are affected by the ioniza- 
tion and recombination processes. Most analytical stud- 
ies in the linear regime adopt the Boussinesq approx- 
imati on where ion arid neutral densities are kept con- 
stant (IBlaes fc Balbu!j[T99llKunz fc Balbusl[200l IDeschI 
I2004D . These studies show that the growth of MRI 
is suppressed when the collision frequency of a neutral 
with the ions falls below the orbital frequency. In the 
mean time, when both vertical and azimuthal field is 
present, unstable modes always exist due to the effect 
of AD, and these unstable modes requir e non-zero ra - 
dial wavenumbe rs (iKunz fc Balbuj [20p: IDeschI [200l . 
IBlaes k Balbuj (|1994[) also studied the effect of ioniza- 
tion and recombination with compressibility (for vertical 
propagating waves), and found that the presence of az- 
imuthal and radial field allows the coupling to acoustic 
and ionization modes, and the azimuthal field tends to 
stabilize the flow when the recombination time is not too 
long compared with dynamical time. 

The effect AD on the MRI in the non-h near regime 
was first studied by iMac Low et al.l (jl995f ). They im- 
plemented and tested AD in the "strong coupling" limit 
(see below) in the ZEUS code and performed simula- 
tions with net vertical flux for various ion-neutral cou- 
pling stre ngths. Their re s ults c onfirmed the stability 
analysis of IBlaes &: Balbu"i (|I994[ ). but their simulations 
are only two-dimensional and did not follow the evolu- 
tion much beyond the li near stage. In another study, 
iBrandenburg et aTl ()I995[ ) included the effect of AD (also 
in the strong coupling limit) in their three-dimensional 
simulations of a local, vertically stratified disk. They 
found that turbulence remains self-sustained in a case 
where AD time is long compared with orbital time, al- 
though reduced in strength, and in another case where 
the AD time was set comparable to fi^^, turbulence de- 
cayed. 

A systematic study o n the non-linear evolution of 
MRI with AD is done bv lHawlev k Stond (|1998f ) (here- 
after HS) using three-dimensional (3D) numerical simula- 
tions. They used the two-fluid approach without consid- 
ering the ionization-recombination processes, therefore 
ions and neutrals obey their own continuity equations. 
Both net-vertical and net-toroidal magnetic configura- 
tions were considered. They found that the system be- 
haves like fully-ionized gas when the ion-neutral collision 
frequency is greater than lOOfi, while ions and neutrals 
behave independently when the collision frequency fall 
below O.OIO. The amplitude of magnetic field at satura- 
tion is proportional to the ion density when it is much 
smaller than the neutral density. The two-fiuid approach 
adopted by HS is valid when the recombination time scale 
is long compared with the dynamical time. However, AD 
in PP Ds is in general in the "strong coupling" limit (jShul 
Il991t) . Two conditions must be satisfied in this limit: 1. 
The ion density pi is negligible compared with the neu- 



tral density p„. 2. The electron recombination time 
must be much smaller than the orbital frequency 51. In 
this limit, the ion density is purely determined by the 
ionization-recombination equilibrium with the neutrals, 
and the two-fiuid formulation is simplified into a single- 
fluid formalism (for the neutrals). In PPDs, condition 
1 is al ways satis fied, and we will show in a companion 
paper ()Baill201ll ) that condition 2 is almost always sat- 
isfied. 

In this paper, we conduct three-dimensional (3D) lo- 
cal shearing-box simulations to explore the effect of AD 
on the non-linear evolution of MRI in the strong cou- 
pling limit. This is conceptually different from the sim- 
ulations performed by HS in that the ion density does 
not obey continuity equation, and is set by the neutral 
density due to chemical equilibrium. Effectively, this al- 
lows the coupling of the MRI with acoustic and ioniza- 
tion modes, which lead s to more complicated interactions 
(jBlaes k, Balbud [19941 ). Moreover, our simulations cor- 
respond to the limit where the ion density is negligibly 
small (i.e., / — >■ in HS), which is difficult for two-fluid 
simulations due to the stiffness of the equations. In the 
strong coupling limit, there is only one controlling pa- 
rameter, namely, the ion- neutral collision frequency ^pi. 
We perform three sets of simulations with net vertical 
flux, net toroidal flux and both. In each group of runs, 
we systematically vary ^pi as well as the strength of the 
net field. Our main goal is to study the conditions under 
which MRI turbulence can be self-sustained or is sup- 
pressed due to AD. In addition, we study the properties 
of the MRI turbulence in the AD dominated regime. 

This paper is structured as follows. In Section [5] we 
provide the formulation of AD in the strong coupling 
limit, describe the numerical method and code test prob- 
lems. A series of numerical simulations on the non-linear 
evolution of MRI with AD are presented and analyzed 
in Section O from which we discuss the condition under 
which MRI turbulence is sustained or suppressed as well 
as the properties of the MRI turbulence in AD domi- 
nated regime. We conclude and bricfiy discuss various 
implications In Section [5l 

2. FORMULATION AND NUMERICAL METHOD 

2.1. Ambipolar Diffusion in Weakly Ionized Plasma 

In weakly ionized plasma, the inertia of the ionized 
species is negligible. The gas dynamics can thus be de- 
scribed by single-fluid equations for the neutrals, modi- 
fled by non-ideal MHD effects. The effect of ambipolar 
diffusion (AD) derives from the relative motion between 
the ions and the neutrals. Since the ion inertia is negligi- 
ble, the ion velocity is determined by the balance between 
the Lorentz force and the ion-neutral collisional drag 

J X B 

= jpnM(v, - v) , (1) 

Hi 

where m;, n.i and Vi are the ion mass, number density and 
velocity respectively, p, v are the density and velocity of 
the neutrals, J , B are the current density and magnetic 
field vectors, 7 = {(7v)/{mn + rrii), with (av) begin the 
rate coefficient for momentum transfer between the ions 
and the neutrals and to„ the neutral mass. The magnetic 
field is effectively carried by the ions, thus the induction 
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equation now reads 

dB 

'dt 



V X (■?; X B) + V X 



{J X B)x B 



An 

V X {v X B) V X 



(2) 



where pi is the ion mass density, va = B / ^Anp is the 
Alfven velocity, Jj. is the component of the current den- 
sity that is perpendicular to the direction of the magnetic 
field. The second term on the right hand side is the AD 
term, from which the ambipolar diffusivity can be defined 
as TjA = v\l^pi. The rest of the fluid equations remain 
the same as the ideal MHD equations. In particular, the 
momentum equation is unchanged since the neutrals ef- 
fectively feel the magnetic tension and pressure by the 
ion-neutral collisions. 

The above derivation is strictly valid when electrons 
and ions are the only charge carriers, and all ions have 
the same mass. Nevertheless, the effect of multiple ion 
and charged grain species can be combined into an ef- 
fective AD coefEcienlQ. AD becomes the dominant non- 
ideal MHD effect when the gyro-frequency of both the 
ions and the electrons are higher than their collision fre- 
quency with the neutrals (i.e., both ions and electrons 
are coupled to the magnetic field). In practice, this cor- 
responds to regions with low density and strong magnetic 
field. 

The effect of AD in rotating disks (with angular 
frequency 17) is characterize d by the parameter Am 
(IChiang fc Murrav-Clavl[200l : 



Am . 



(3) 



which is the number of effective collisions for a neu- 
tral molecule/atom with the ions in a dynamical time 
l/f2. Physically, we see that Am = v\/r]A^, which is 
equivalent to the Elsasser number for Ohmic resistivity 
A = v\/rj^. Therefore, Am measures the ratio of the 
AD time scale over the critical wavelength va/^ and the 
dynamical time scale. 

In this paper, we study AD in the strong coupling 
limit ([Shu 199ll): the electron (ion) recombination time 
is much shorter than the dynamical time so that 
the ion density pi is determined by the local thcrmo- 
dynamical quantities (density and temperature). The 
strong coupling limit is widely applicable in protop lan- 
etary disks, and we show in our companion paper (jBail 
120 lit) that the recombination time is typically at least 
one order of magnitude smaller than the dynamical time, 
even in the disk coronal regions. Using the strong cou- 
pling limit, we assume the ion density depends on the 
neutral density p in the form of 



, P 
P« — 



(4) 



where p^o Eind po are the reference density of the ions 
and the neutrals. A simple-minded calculation of the 

^ A more generalized derivation of all non-ide al MHD ter ms in 
the ind uction e quation is given and discussed in IBail I I2011I ) (and 
see also lWardlel 1 1199 9. .2007i) ), which also includes Ohmic resistivity 
and the Hall term. 



ionization-recombination equilibrium gives = Cpf, 
which yields = 1/2. In the equation ^ is the ioniza- 
tion rate, of the order 10~^^ s~^, C is the effective rate 
coefficient for recombination, of the order 10^-^ cm^ s~^ 
g~^ in the absence of grains (jBlaes fc Balbusl[T99l . In 
reality, the recombination process is complicated by a 
complex network of gas phase and grain phase chemical 
reactions, and we can address the complications by ex- 
ploring different values of v. Being a two-body process 
in general, and one expects < < 10 In fact, our sim- 
ulations show that the properties of the MRI turbulence 
is insensitive to v (see Section fS.l.Sp . 

2.2. Numerical Method 

We use Athena, a higher-order Godunov MHD 
code with constrained transport technique to enforce 
the divergence-free constrain t on the magnet ic field 
(IGardiner fc Stonel[2(M [20081 : IStone et al.l[2008h for all 
calculations presented in this pa per. Non-ideal M HD 
terms including Ohmic resistivity (jDavis et al.ll20ldl ). the 
Hall term (in progress) and AD (this paper) have been 
developed for Athena. We consider a local patch of a 
Kepl erian disk using the standard s hearing-box formal- 
ism (jGoldreich &: Lvnden-Belilll965[ ). which adopts a lo- 
cal reference frame at a fiducial radius corotating with 
the disk at orbital frequency f2. In this frame, we write 
the MHD equations with AD in a Cartesian coordinate 
system as 



dp „ 



{pv) ^ , 



^ + v.(p.-. 



dB 

'dt 



V X 



T)=p 



V X B 



2v xn + sn'^xx 



{J xB)xB 



where T is the total stress tensor 
T = (P + BVStt) I - 



CJiPiP 



B^B 



47r 



(5) 



(6) 



(7) 



(8) 



I is the identity tensor, P is the gas pressure, x, y, z are 
unit vectors pointing to the radial, azimuthal and vertical 
directions respectively, where $7 is along the z direction. 
Disk vertical gravity is ignored and our simulations are 
vertically unstratified. We use an isothermal equation of 
state P = pCj, where Cs is the isothermal sound speed. 
Periodic boundary conditions are used in the azimuthal 
and vertical directions, while the radial boundary condi- 
tions are shearing periodic as usual. 

An orb ital advection scheme (jMassetl 120001 : 
iJohnson et al. 2008 i) has been implemented in Athena 
(jStone fc Gardineii l2010f) . It splits the dynamical 
equations into two systems, one of which corresponds to 
linear advection operator with background flow velocity 
3nx/2y and the other evolves only velocity fluctuations. 
The orbital advection scheme not only accelerates the 
calculation in large box size by admitting larger time 
steps, but also makes the calculation more accurate. 



For example, Figure la and Fig ure 3 of lBai fc GoodmanI if^OOgP 
illustrate the dependence of electron abundance oc pi on the gas 
density in various chemistry models v^ith and without dust grains, 
and < < 1 holds in essentially all circumstances. 
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The last term on the right hand side of equation ([7]) 
represents the AD term, where pi is the effective ion 
density, approximated by equation (|4]). Non-ideal MHD 
terms (e.g., the AD term) arc implemented in Athena 
in an operator-split way. Over one time step, one first 
solves the induction equation using only the non-ideal 
MHD terms, and then solves ideal MHD equations. The 
induction equation with AD term is a diffusion equation 
and is evolved by a fully explicit forward-Euler method. 
The divergence free condition is maintained using con- 
strained transport: we calculate the AD electromotive 
force defined at cell edges by interpolating the magnetic 
field and current density to such locations. The method 
is stable, with time step constrained to be proportional 
to grid size squared. The overall method is first order 
accurate in time. 

Numerical method including the AD term has been fre- 
quently implemented in the framew ork of two- fluid mod- 
els (e. g.-l^thllTQQllStOT ic 1997; Smit h fc Mac Low 1997; 
FalS 20031; iLi et al.l 12006,; ,0'Sullivan fc Downcs, ,200a 



20071 iTillev fc Balsara]l20'08h . with the main applications 
on star formation, turbulence and shocks in the inter- 
stellar medium. However, single-fluid models for AD 
in the strong coupling limit relevant to weakly ion ized 
disks ()Brandenburg et al.lll995HMac Low et al.lll995[) are 
relativ ely less studied numericallv. Recentlv. iChoi et al.l 
()2009D described an explicit scheme for incorporating AD 
in the strong coupling limit in an MHD code that is sec- 
ond order accurate, and use super time-stepping to ac- 
celerate the calculation when the AD coefficient is large. 
As we find MRI is suppressed at moderately large AD 
coefficients (see Section O , super time-stepping is not 
needed for the purpose of this paper. 

2.3. Code Tests 

In this subsection, we conduct two test problems to 
examine the performance of our implementation of the 
AD term. For all these tests, we consider non-rotating 
system and turn off the shearing box source terms on the 
right hand side of equation ([6]). 

2.3.1. Isothermal C-type Shock Test 

The effect of ambipolar di ffusion is best manifested in 
C-type shocks ()Drainelll980[ ). which is a shock with con- 
tinuous transitions consequent of the AD. For the pur- 
pose of the c ode test, here we cons ider the isothermal 
C-type test bv lMac Low et abl ()1995D . which has become 
a standard test problem for AD. 

We consider steady-state (d/dt — 0) shock and work in 
the shock frame, with upstream gas density po moving at 
velocity Vg. The upstream gas is threaded by a uniform 
magnetic field Bq that lies at an angle 6 to the velocity. 
Let Vs be in the x direction, and Bq be in the x — y plane. 
For a continuous shock, the jump conditions reduce to 
d/dx = 0. Assuming the ion density pi is constant {v = 
0), the equations that describe the C-type shock read 



PVxVy 



PVx 

B 



PoVs 



-lL^p,iv'^+cl) + ^ 



BxBy 

An 



BxBoy 

in 



(9) 
(10) 

(11) 



VxB„ 



VyBx 



B^ dBy 

Anp^Pi dx 



v.Bt 



Oy 



(12) 



Note that Bx = Bqx is constant. 

The shock is characterized by three dimensionless pa- 
rameters: the sonic Mach number M = Vg/cs, the Alfven 
Mach number A = Vs/va (where v\ = Bg/inpo), and 
the angle 9 of the magnetic field with the upstream 
flow. The characteristic length scale of the problem is 
given by L = va/jpi- We further define D = p/po, 
and b = By/Bo. After some algebra, we arrive at 
a dimensionless first order differential equation for D 
pvlac Low et aLl ll995D 



1 



1 



dD 



AP J d{x/L) 



A 



b-D 



sm ( 



A^ 



■ cos 



(13) 



One can numerically integrate this ordinary differential 
equation to obtain the C-type shock profile. In Figure [T] 
we show a semi-analytical solution for M = 50, A = 10 
and = n/A obtained by using a 4th order Runge-Kutta 
method. 

To use this solution as a code test, the shock is set 
to be aligned with the grid in the x direction. We use 
outflow boundary conditions in this direction. In multi- 
dimensional tests, periodic boundary conditions are used 
in other directions. The shock solution should be station- 
ary (i.e., a standing shock), thus we evolve the solution 
for sufficiently long time (^ bL/cs) and compare to the 
initial conditions. In Figure [U we further show the ab- 
solute error of the shock profile compared with the semi- 
analytic solution. Since the shock is grid-aligned, ID, 2D 
and 3D tests essentially produce the same result. In our 
tests, the grid resolution is chosen to be 2 and 4 cells 
per L0 We see that our code very accurately resolves 
the structure of the C-type shock using only a few cells 
per L. The main source of the error lie in the region 
where density and velocity profile vary quickly. In our 
comparison, the position of the shock is fixed at the ini- 
tial place, while in reality, the shock position can shift 
slightly during numerical relaxation. 

2.3.2. Damping of MHD Waves 

Linear MHD waves are damped due to AD. Because 
the exact eigenvectors in the AD regime can be very 
complicated, we initialize the problem with ideal MHD 
wave eigenvectors and measure the damping rate. This 
means the initial conditions are a linear superposition 
of more than one eigen-mode, but the averaged damp- 
ing rate should approach the analytical value for a single 
mode as long as the AD coefficient is sufficiently small. 

The analytical damping rate for various MH D waves 
due to AD can be found in and/or derived from lBalsaral 
(|1996f l. as we summarize below. The damping rate of the 
Alfven wave is given by the solution of 



= k^vi cos 



1 - i — 



(14) 



^ In comparison, IMac Low et al.l l l 19951) ac hieved comparable 
accuracy as ours using 5 and 10 cells per L, IChoi et al.l ||2009|'| 
achieved similar or better accuracy at 6.4 cells per L. 
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Fig. 1. — The profile of a C-typo shock with M = 50, A = 10 and 8 = n/A. The upper panels show the semi-analytical solution of 
gas density p, perpendicular velocity Vy and perpendicular magnetic field By. p and By are normalized to their upstream values, and Vy 
is normalized to the sound speed. The lower panels show the corresponding absolute errors (same units as the upper panels) from our 
numerical simulations, with resolutions of 4 cells per L (black solid) and 2 cells per L (blue dash-dotted). 



where = 7Pi, va is tlic Alfvcn velocity, 9 is the angle 
between the magnetic field and the wave vector. The 
damping rate of fast and slow waves can be obtained by 
solving the quadratic equation 

(^2 _ ^2)(^2 _ „2) ^ i(^2 „ k^Dk^vl- = , (15) 

where Vf and Vg arc the fast and slow magnetosonic 
speeds. 

When cj <C the damping rate is small and can be 
found by expanding the ideal MHD dispersion relation 
to powers of uj/uja, and to the first order, we find for the 
damping rate of the Alfvven wave 

^^^levlcosH _ ^^^^ 

2 LOa 

The damping rate for fast and slow magnetosonic waves 
are 

We perform the linear wave damping test in ID, 2D 
and 3D. In ID, the wave is grid-aligned, whose wave 
length equals 1 in code unit. In 2D and 3D test problems, 
the wave vectors arc not grid-aligned, with box sizes cho- 
sen such that the wave length is also 1 [in 2D, the box 

size is {Vb, ■\/5/2) and in 3D, the box size is (3, 1.5, 1.5)]. 
We use isothermal equation of state with = 1. The 
background gas density is po = 1- As before, we choose 
v = 0. In ID, the wave vector is along the x direction, 
and the adopted magnetic field is Bq^ = 1.0, Bgy = \f2 
and Bqz = 0.5. In 2D and 3D the background magnetic 
field vector is rotated with the wave vector accordingly 
while keeping 9 the same as in ID (and a vector poten- 
tial is used to initialize the wave in order to preserve 
the divergence free condition). From the above we get 
VA = \/13/4, Vf = 2 and Vs = 1/2, therefore we find 
= 27r^QA.. T/ = 5.2Tr^QA and = I.Stt^Q^. 



In practice, we adopt uja = 100 and run our simulation 
to i = 5. By default, the grid size is 32 in ID, 64 x 32 
in 2D, and 64 x 32 x 32 in 3D. Accounting for the box 
size in each dimension, the effective resolution, charac- 
terized by number of cells per wavelength, is 32, 28.6 
and 21.3 in ID. 2D and 3D respectively. The results are 
shown in Figure [2] From left to right, we show the damp- 
ing curve from ID, 2D and 3D simulations in the solid 
lines, where black, blue and red lines label Alfven, fast 
and slow MHD waves. Dashed lines show the theoretical 
damping curve. We see that the numerical damping rate 
matches very well with the theoretical damping rate. In 
the 3D runs, the damping rate for the Alfven wave is 
slightly faster than expected, but this may be because 
the effective resolution is less. 

Wc have also run the simulations with double and half 
the resolution. With double resolution, the numerical 
damping curves in ID, 2D and 3D cases almost match ex- 
actly the analytical damping curves (besides some small 
oscillations due to the initial conditions). At half resolu- 
tion, however, the numerical damping rate deviates sub- 
stantially (about 15% to 30% at t ~ 5). These results 
indicate that at least 20 cells per wavelength is need to 
accurately capture the effect of AD. 

3. SIMULATIONS AND RESULTS 

In this section we describe three groups of simulations 
and study the effect of AD on the non-linear evolution 
of the MRI. All our simulations are vertically unstrati- 
fied by ignoring the disk vertical gravity with fixed box 
height to be one disk scale height H = Cs/il. We ini- 
tialize our simulations with Keplerian velocity and seed 
density perturbations of 2.5% of the background density 
Po- We consider three different magnetic field geometries 
(net vertical fiux, net toroidal flux, both vertical and 
toroidal flux) as described in the following three sub- 
sections. Since all our simulations contain net magnetic 
flux, they arc not subject to the issue of convergence 
found by (Froman g fc Papaloizoul [20071 ) in zero net-flux 
simulations, and numerical convergence is confirmed in 
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Fig. 2. — The damping of linear MHD waves by ambipolar diffusion. Three panels from left to right show the test results for ID, 2D and 
3D, where in the latter two cases the waves arc not grid-aligned. In each panel, black, blue and red curves show the damping of Alfvvcn, 
fast and slow waves respectively. Solid linos arc the measured damping curve, while dashed lines are the expected damping curve. 



our test simulations (and see Section 13.2.21 for the case 
of net toroidal flux). For relatively small AD coefficient 
(large Am), MRI grows quickly from the seed pertur- 
bations and saturates into turbulence; when the effect 
of AD is strong, however, MRI does not grow from our 
seed perturbations. In such cases, we initialize the sim- 
ulations from a turbulent state which is obtained from 
simulations with relatively large Am (sec individual sub- 
sections for details). 

The most important diagnostics are the volume aver- 
aged (normalized) Reynolds stress, defined as 



(18) 



where the over bar indicates volume averaging, v'y is 
the azimuthal velocity with the Keplerian velocity sub- 
tracted, and the volume averaged (normalized) Maxwell 
stress, defined as 



«Max — 



47rpoCs 



(19) 



The total stress, na mely, the a parameter 
(jShakura fc Sunvaevlll973l ) is a = + Q^Max- We also 
monitor the kinetic and magnetic energy density, which 
characterize the strength of the MRI turbulence. 

The main purpose of this study is to identify the cri- 
terion under which sustained turbulence generated by 
the MRI can be maintained. However, the term "sus- 
tained turbulence" is a somewhat ambiguous concept. 
In the context of 3D shearing box simulations, small- 
amplitude oscillations left from decayed hy drodynami- 
cal turbulence is present (jShen et al.l 120061) . Such os- 
ciUations produce Reynolds stress on the order of 10^^ 
with kinetic energy density on the order of 10~^, both 
in normalized unit poc^- Such oscillations are likely to 
associate with linear modes in the shearing sheet with 
the origin of acoustic an d/or nearly incompressible iner- 
tia waves (jBalbusI [20031 ). Being a conservative Godunov 
MHD code, the Athena code preserves the amplitude of 
these waves without much damping. Therefore, through- 
out this paper, the level of turbulence we are interested 



in are those whose time and volume averaged kinetic en- 
ergy density Ek = {pv"^) is on the order of 10~^pocl or 
higher, and/or whose total stress a is no less than 10^**. 
Meanwhile, analysis of all our simulations show that the 
threshold where the MRI turbulence can be marginally 
self-sustained is roughly at the same level, (see Section 
13.2.41 for further discussion). 

Our simulations are run for at least 24 orbits (150f2^^). 
A period of 24 orbits is sufficiently long for the MRI to 
saturate from initial growth, which typically occurs in 
5 — 10 orbits, or for the restart runs to reach a steady 
state, which typically occurs in 10 — 15 orbits. Our time 
averaged quantities are mostly taken from after about 16 
orbits (since 10017"^) unless otherwise noted. Although 
a time average over 8 orbits (50r2~^) is relatively short, 
it is sufficient for our purpose to judge whether MRI tur- 
bulence can be self-sustainccfl. Many of our simulations 
are run for 48 orbits or longer where better statistics on 
the turbulence properties can be obtained. 

3.1. Net Vertical Flux Simulations 

In the first group of simulations, we choose the ini- 
tial field configuration to be uniform along the vertical 
axis z, characterized by the plasma /3o = ^Pq/Bq, where 
Pq = poCg is the background pressure and Bq is the initial 
field strength. The vertical flux is conserved numerically 
by remapping the toroidal component of the magnetic 
field in th e ghost zones of th e radia l boundaries (see Sec- 
tion 4 of IStone fc Gardined (I2010D for details). For aU 
simulations, we fix the box size to be 4iJ x AH x H 
in the radial, azimuthal and vertical dimensions, with 
fixed grid resolution at 64 cells per H. We have cho- 
sen a relatively large radial b ox size (4iJ), as suggested 
by iPessah fc GoodmanI ()2009D , which is needed to prop- 



^ IWinters et al.l II2003I ) found that more than a few hundred or- 
bits as are required to accurately measure the properties of the 
MRI turbulence in ideal MHD. This conclusion is based simula- 
tions with radial box size being H, while the radial box size in 
about half of our simulations is AH, which reduces the time fluc- 
tuations. Also, our Figures l3l and [8l show that the fluctuations in 
the Maxwell stress are less severe in the presence of AD than in 
the ideal MHD case. 



Effect of Ambipolar Diffusion on MRI 



7 



erly capture the parasitic modes to break the channel 
mode into turbulence. It also help substant ially reduce 
the intermittence of the MRI turbulence (|Bodo et alj 
l2008f l. We note that for local unstratified net verti- 
cal flux MRI simulations without explicit dissipation, 
turbulence propertie s converge at about 32 cells per H 
(|Hawlev et al.lfl995[ ). The grid resolution in our sim- 
ulations is two times higher, thus we expect numerical 
convergence. 

All our net vertical flux simulations arc listed in Table 
[TJ We first perform a fiducial set of simulations with fixed 
/3o = 400. We choose a series of Am values, ranging from 
1000 down to 0.1, and study the critical value of Am 
below which MRI turbulence is no longer self-sustained 
(Section l3.1.ip . In the next, we vary the net vertical flux 
by setting = 100, 1600 and lO" and run a number 
of simulations around Am = 1 to study how the critical 
value of Am is affected by the vertical flux fScction l3.1.2|) . 
Moreover, in Section [3.1.3l we briefly investigate the effect 
of v by varying v from the fiducial value 0.5 to (run Z5a) 
and 1 (run Z5b) (see equation ([3])). Finally, we discuss 
the properties of the MRI turbulence in the presence of 
AD (Section [3X31). 

Our choices of the net vertical flux derive from the 
linear dispersion relation of the MRI as well as phys- 
ical considerations. In the case of ideal MHD, the 
wavelength for the fastest growing linear MRI mode is 

given by X/H = QASP^^'"^ (IHawlev et al.l[l995l) . For 
= 100, 400 and 1600, our vertical box size of H fits 
1, 2 and 4 most unstable wavelengths respectively in ideal 
MHD. The ideal MHD dispersion relation is considerably 
modified when Am ^10. Unstable modes exist for wave - 
length longer than the critical wavelength (|Wardlelll999() 

The wavelength for the most unstable mode Am is about 
twice larger. An approximate fitting formula that is ac- 
curate within 2% for all values of Am is 

^.10,26(l + J-, + -^-0,2e)"V"',(21) 

where e = Am/{\ + Am). Note that for pure verti- 
cal magnetic field and vertical wavenumber, the linear 
dispersion relation for Ohmic resis tivity is exactly the 
same as that for AD () Wardld [l999l ) . with Am replaced 
by the Elsasser number A = v\^/rj^. For /3o = 400, the 
most unstable wavelengths at Am = 1,1/3 and 0.1 are 
A = 0.87H, 1.72H and 5.18H respectively. Clearly, the 
most unstable mode docs not fit into our simulation box 
when Am = 0.33, and no unstable modes fit into the 

— 1/2 

box for Am = 0.1. Since A oc /3q , these modes do fit 
into our simulation box as we increase /3o to 1600 and 
10^ respectively. In the mean time, since AD tends to be 
important in the more stron gly magnetized upper lay ers 
of the protoplanetary disks (jWardld [2007L lBajfl201l( i. it 
is also interesting to study whether the MRI turbulence 
can be sustained when /3o is relatively small, even if the 
most (or all) unstable modes do not fit into our simula- 
tion box. We have not run simulations with a taller box 
since we do not include vertical stratification. 



TABLE 1 
Net vertical flux simulations. 



Run 


Am 


/3o 


V 


Orbits 


Restart^ 


Turbulei 


Zl 


1000 


400 


0.5 


48 


No 


Yes 


Z2 


100 


400 


0.5 


24 


No 


Yes 


Z3 


10 


400 


0.5 


24 


No 


Yes 


Z4 


3.33 


400 


0.5 


24 


No 


Yes 


Z5 


1 


400 


0.5 


43 


No 


Yes 


Z6 


0.33 


400 


0.5 


24 


Z5 


No 


Z7 


0.1 


400 


0.5 


24 


Z5 


No 


Z3s 


10 


100 


0.5 


24 


No 


Yes 


Z5s 


1 


100 


0.5 


24 


No 


Yes 


Z6s 


0.33 


100 


0.5 


24 


Z5s 


No 


Z7s 


0.1 


100 


0.5 


24 


Z5s 


No 


Z3w 


10 


1600 


0.5 


24 


No 


Yes 


Z5w 


1 


1600 


0.5 


24 


No 


Yes 


Z6w 


0.33 


1600 


0.5 


48 


No 


Yes 


Z7w 


0.1 


1600 


0.5 


24 


Z5w 


No 


Z6e 


0.33 


10* 


0.5 


48 


No 


Yes 


Z7e 


0.1 


10* 


0.5 


48 


Z6c 


Yes 


Z5a 


1 


400 


0.0 


24 


No 


Yes 


Z5b 


1 


400 


1.0 


24 


No 


Yes 
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Box size is fixed at 4_ff X 4_ff X H, grid resolution is 64 cells per H. 
^ whether simulation is initiated by restarting from a turbulent 
run. 

whether turbulence can be self-sustained. 

3.1.1. A Fiducial Set of Runs 

As the fiducial set of runs, we fix /Sq = 400, and run 
7 simulations with different Am values (see Table [T]). la- 
beled from Zl with Am ~ 1000, which essentially cor- 
responds to the ideal MHD case, to Z7 with Am = 0.1, 
where the evolution of magnetic field is dominated by 
AD. Our scan of Am is more narrowly sampled near 
Am = 1 where the transition is expected to occur. In 
Figure [3l we show the time evolution of the Maxwell 
stress from the fiducial set of runs. We find that for 
Am > 1, the growth of the MRI from linear perturba- 
tions leads to vigorous MRI turbulence, while for the two 
runs with Am < 1, MRI either does not grow from the 
initial vertical field (Z7), or grows too slowly (Z6), since 
the most unstable modes do not fit into our simulation 
box. Therefore, for these two models, we start the simu- 
lations from the end of run Z5 (Am = 1), which is turbu- 
lent, and reset Am to be 0.33 and 0.1 respectively. Nev- 
ertheless, turbulence continues to decay throughout the 
span of our simulation in run Z7. Run Z6 is a marginal 
case where turbulence is neither fully sustained nor de- 
cayed continuously (see discussions below). 

We first look at run Zl to Z5 with Am > 1. The ini- 
tial growt h of the MRI is due to t he axisymmetric chan - 
nel mode (|Goodman fc Xul [1991 Fessah fc Chai^lMl . 
The mode becomes non-linear (producing an overshoot 
in the Maxwell stress up to 1 in the ideal MHD case) 
until broken down by secondary parasitic modes to pro- 
duce turbulence. In the turbulent state, it is evident 
that the Maxwell stress monotonically decreases as Am 
decreases, analogou s to the Ohmic case (ISano et al.lll998l : 
[Fleming et al|[2000[ ). In Table [5] we list the general prop- 
erties of the turbulence from all our vertical net flux sim- 
ulations. The quantities are averaged over space and 
time after saturation (after 10017"^). The total stress 
a w 0.2 in run Zl, which agrees with the ideal MHD 
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Fig. 3. — The evolution of Maxwell stress in our fiducial set of net 
vertical flux runs (from top to bottom: Zl, Z2, Z7) normalized 
to CsH . For models Z6 and Z7, the simulations are initiated from 
the end of run Z5. 



case (|Hawlev et al.l[l995l ). It drops slowly with decreas- 
ing Am when Am 3> 1, but very rapidly when Am is 
around 1. Moreover, as Am decreases, the ratio of kinetic 
to the fluctuating part (with background field Bzo sub- 
tracted) of the magnetic energy increases (see also Fig- 
ure [S]). Similarly, the ratio of Reynolds stress to Maxwell 
stress increases. 

As we have discussed before, the most unstable mode 
does not fit into our simulation box for runs Z6 and Z7, 
and our simulations initiated from a turbulent state also 
show no sign of sustained MRI turbulence. This is not 
surprising for run Z7, where no unstable MRI mode even 
exists in the simulation box. Our run Z6 is a marginal 
case, where a wavelength of H is only slightly larger than 
the critical wavelength for instability (Ac — 0.89i?) but 
far from the most unstable wavelength (Am = 1.7277). 
This explains the long-term variations in the Figure [3] 
since the growth rate is only slightly larger than zero. 
Our analysis in Section 13.2.41 indicates that although 
some non-zero stress close to a = 10"'' is maintained 
in the simulation, it is unlikely to be due to the MRI 
turbulence. In real disks, one may expect sustained tur- 
bulence to be supported at Am sa 0.3 if it were at the 
disk midplane, where the density variation over one H 
above and below the midplane is not significant. In the 
upper layers, /3o may fall off substantially over one H , 
which strongly stabilizes the flow. In sum, we see from 
our fiducial set of net vertical fiux simulations that pres- 
ence of turbulence mainly depends on whether the most 
unstable mode of the MRI fits into the simulation box. 
This aspect will be further explored in the next subsec- 
tion [3X1 

The results reported above are qualitatively different 
from those observed in HS using two-fluid simulations. 
One may compare our results with Table 3 of HS, where 
the Am value for their four runs Z24, Z17, Z25, Z28 are 
0.11, 1.1, 11.1 and 111 respectively. The total stress a in 
these simulations does not scale monotonically with Am, 
and in particular a is on the order of 10~^ when Am is as 
small as 0.11, while in our simulations MRI is suppressed. 



This reflects the difference in the physical assumptions 
about the two approaches. In the two-fiuid limit, ions 
are coupled to the neutrals only via collisions. When the 
ion- neutral collisions are infrequent [Am < 1), the ions 
and neutrals behave as independent fluid: vigorous MRI 
is generated in the ion fluid while the neutrals remain 
quiescent, and the overall a is proportional to the ion 
fraction f = pi/p. In the strong-coupling limit, the ions 
and the neutrals are coupled not only by collisions, but 
also via the ionization-recombination reactions. Since 
the strong-coupling limit requires the recombination time 
to be much smaller than the orbital time, this means 
that ions are continuously created and destroyed on a 
time scale that is much shorter than the time scale for 
MRI to grow. It is this additional chemical coupling 
that quenches the MRI in the ion fluid, and suppresses 
angular momentum transport for Am ^0.1. 

3.1.2. The Effect of Vertical Field Strength 

We select a number of models from the fiducial series 
and rerun the simulations with three additional initial /3o 
values: /3o = 100, 1600 and 10"* (e.g., with magnetic field 
strength two times and half of that in the fiducial mod- 
els, as well as one case with an very weak field). These 
simulations are labeled with an additional letter "s" (for 
strong), "w" (for weak) and "e" (for extremely weak) 
in Table [TJ For strong field simulations with /3o = 100, 
the wavelength for the fastest growing mode exceeds the 
vertical box size when Am < 10. When Am < 1, there 
are essentially no unstable mode in the simulation box. 
Runs Z3s and Z5s are initiated from seed perturbations, 
while runs Z6s and Z7s are initiated from the turbulent 
state at the end of run Z5s. For weak field runs with 
(3q ~ 1600, on the other hand, the most unstable mode 
can be fitted into the simulation box for all runs with 
Am > 0.3. No unstable mode is fitted into the simu- 
lation box when Am = 0.1, and as before, Run Z7w is 
initiated from turbulent state from Z5w to test whether 
turbulence can be sustained. Our /3o = 10^ simulations 
allows the most unstable wavelength to be fitted in the 
simulation box at small Am. We conduct two runs in 
this case with Am = 0.33 (Z6e) and Am = 0.1 (Z7e). 
Run Z7e is initialized from the turbulent state in Z6e to 
avoid the extremely long time in the linear growth stage. 
Time averaging in runs Z6w, Z6e and Z7e are taken since 
t = 200n~^ (time averaging in other runs are taken since 
t = 10017-1 by default). 

We find from our simulations that sustained MRI tur- 
bulence is present in all models except Z6s, Z7s and Z7w. 
In particular, the MRI turbulence can be self-sustained 
even the Am value is as small as 0.1, provided that the 
net vertical field is sufficiently weak. These results con- 
firm our speculation in the fiducial set of simulations that 
the MRI turbulence is self-sustained as long as unstable 
MRI modes fit into the simulation box. 

The diagnostic quantities from time and volume aver- 
aged quantities in the turbulent state from the weak and 
strong field series of runs are also listed in Table[2j We see 
that for Am = 10, the averaged kinetic energy, Reynolds 
and Maxwell stress monotonically decreases with increas- 
ing /So- Although not all our simulations are not run 
long enough for these quantities to be measured accu- 
rately, the trend is significant enough and indicates that 
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Fig. 4. — The time and volume averaged total stress a. from 
all our net vertical flux simulations that sustain MRI turbulence. 
Simulations with different net vertical flux, characterized by the 
plasma /3o, are labeled by different symbols and colors. The ar- 
rows above the symbols indicate (for each /3o as represented by the 
symbol) the range of Am where the most unstable wavelength is 
smaller than H. The dashed line connecting the symbols repre- 
sents the maximum value of stress attainable from net vertical flux 
simulations. 

the MRI saturate at a higlier level with higher net ver- 
tical flux (small /3n) , in a greement with the ideal MHD 
case (jHawlev et al.lfl995() . For Am = 1, the monotonic- 
ity trend is still present by comparing our fiducial run 
Z5 and the weak field run Z5w. The saturation level of 
the MRI turbulence in the strong field run Z5s is weaker 
than that for run Z5. This is most likely because the 
most unstable mode does not fit into our simulation box 
(but some less unstable modes fit) in run Z5s. The mono- 
tonicity trend further preserves at Am = 0.33, where the 
kinetic energy density and total stress from run Z6w is 
larger than those from run Z6c by about a factor of 2. 

We summarize the main results from the net vertical 
flux simulations in Figure H) Shown arc the total stress 
a from all the simulations that the most unstable mode 
is properly resolved so that the MRI turbulence is self- 
sustained and a reliable value of a can be obtained. As 
discussed before, at a fixed /3o, there exists a critical value 
of Am below which the most unstable wavelength would 
exceed H and the mode tend to be suppressed due to 
the vertical stratification. This is effect is illustrated by 
the colored arrows in the Figure. Equivalcntly, for turbu- 
lence to be sustained at small Am, /3o must be sufficiently 
large /3o ^ as can be obtained from equation 

(|2T|) . At a given Am, since the stress a monotonically 
increases with the net vertical flux, there exist a max- 
imum stress, corresponding to the largest allowed net 
flux (smallest allowed value of /3o). This maximum value 
of a as a function of Am is illustrated in the dashed 
line by connecting results from runs Z7e, Z6w, Z5 and 
Z3s. We see that the maximum a drops by a factor of 
about 40 from the ideal MHD case to Am = 1, and an- 
other factor of about 60 as Am decreases to 0.1. By 
extrapolating this trend, we expect the MRI turbulence 
can be self-sustained for arbitrarily small value of Am, 



as long as the background magnetic field is sufficiently 
weak. Nevertheless, the turbulence would seem to be 
too weak (a < 10~^) to produce significant amount of 
angular momentum transport as required by most astro- 
physical disks. 

3.1.3. The Effect of u 

The parameter v reflects the sensitivity of how the 
AD coefficient depends on gas density (see Equation 
((4])). Most of our simulations arc run with fixed value of 
V — 0.5, while v can in principal span a range from to 
1. The significance about the effect of v largely depends 
on the level of density ffuctuation in the MRI turbulence. 
In Table [21 wc list the rms density fluctuation relative to 
the background gas density from all vertical net flux runs 
(see column {5p) / po)- The rms density fluctuation in the 
ideal MHD case (run Zl) is relatively large, up to 0.3, 
and the largest and smallest densities reach about 0.2 
and 4 times the background density. Since AD reduces 
the saturation level of the MRI turbulence, the density 
fluctuations become smaller as Am decreases. This fact 
undermines the importance of v: when the effect of v 
may be important (large density fluctuations), AD only 
plays an insignificant role in the MRI turbulence (large 
Am); when AD strongly affect the MRI turbulence (small 
Am), the density fluctuation becomes much smaller and 
v is much less likely to be important. This above implies 
that variations in the value of v should not have a major 
impact, and in particular, the critical value of Am be- 
low which MRI is suppressed is unlikely to be altered by 
different choices of i'. 

To confirm our expectations, we perform two addi- 
tional runs with the same initial conditions as run Z5 
(Am = 1, /3o = 400), but set to be and 1 respectively. 
These two runs are named Z5a and Z5b. We see from Ta- 
ble [2] that the turbulence properties from these two runs 
are essentially identical to those in run Z5. Even though 
our time averages arc taken over relatively short periods, 
the deviations are generally within 10%. This is under- 
standable since the density fluctuations in these runs are 
as small as 0.07. It appears certain that the value of ly 
only plays a very minor role in the MRI turbulence in 
the strong coupling limit. 

3.1.4. Properties of the MRI turbulence with AD 

Besides the general properties of the MRI turbulence 
listed in Table [5J we study two other aspects of the MRI 
turbulence with AD. 

First, we study the power spectrum density (PSD) of 
magnetic and kinetic energies by Fourier analysis. The 
Fourier analysis in the shearing periodic system is per- 
formed by the remapping technique before and after 
Fourier trans f ormat ion, as described in Section 2.4 of 
iHawlev et al.l ()1995D . Although the PSD is anisotropic 
in fc— space, it would be bencflcial to plot the PSD in 
one di mensional form by s ome averaging procedure. Fol- 
lowing [DwisIer^I] (j2010f ). we compute shell- integrated 
power spectrum of the magnetic held = 47rfc^|JB(fc)p, 
where |-B(fc)p denotes the average of |S(A;)p over shells 

of constant k = \k\, and B{k) = / B{x)e~^''''"d^x/V is 
the Fourier transform of B{x). Here V is the volume 
of the simulation box. The Fourier transformation is of 
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TABLE 2 

Time and volume averaged quantities in net vertical flux simulations. 



Run 








Ek.y 








Ek 




Em,x 






Em.v 






Em.z 




{opj/po 




a-Rc 




OMax 




a 




Zl 


7.0 X 


10- 


s 




0.10 




2.9 X 


10- 




0.20 




8.4 


X 


10- 






0.20 




3 


4 


X 


10- 


i 


0.36 


5 


4 X 


10- 




0.17 




0.22 




Z2 


4.4 X 


10" 


2 


6 


5 X 


10- 


2 


1.9 X 


IQ- 


2 


0.13 




4.5 


X 


IQ- 


2 


8 


7 X 


10- 


2 


2 


1 


X 


10- 


-2 


0.31 


3 


3 X 


10- 


2 


7.4 X 


10- 


2 


0.11 




Z3 


2.8 X 


10- 


2 


2 


7 X 


10- 


2 


1.4 X 


10- 


2 


6.9 X 


10- 


2 


1.4 


X 


10- 


2 


3 


3 X 


10- 


2 


9 


9 


X 


10- 


-3 


0.20 


1 


6 X 


10- 


2 


2.9 X 


10- 


2 


4.5 X 10- 


2 


Z4 


1.8 X 


10" 


2 


1 


3 X 


IQ- 


2 


8.6 X 


10- 


3 


4.0 X 


10- 


2 


6.5 


X 


10- 


3 


1 


8 X 


10- 


2 


6 


3 


X 


10- 


-3 


0.14 


9 


X 


10- 


3 


1.5 X 


10- 


2 


2.4 X 10- 


2 


Z5 


5.4 X 


10" 


3 


1 


9 X 


IQ- 


3 


3.1 X 


10- 


3 


1.0 X 


10- 


2 


1.2 


X 


10- 


3 


5 


1 X 


10- 


3 


3 


3 


X 


10- 


-3 


0.070 


2 


4 X 


10- 


3 


4.1 X 


10- 


3 


6.5 X 10- 


3 


Z3s 


3.5 X 


10" 


2 


3 


3 X 


IQ- 


2 


9.6 X 


10- 


3 


7.8 X 


10- 


2 


3.4 


X 


10- 


2 


5 


X 


10- 


2 


2 


7 


X 


10- 


-2 


0.18 


2 


9 X 


10- 


2 


4.0 X 


10- 


2 


6.9 X 10- 


2 


Z5s 


3.7 X 


10- 


3 


1 


8 X 


10- 


3 


3.5 X 


10- 


3 


9.0 X 


10- 


3 


1.2 


X 


10- 


3 


1 


2 X 


10- 


3 


1 





X 


10- 


-2 


0.063 


2 


4 X 


10- 


3 


2.2 X 


10- 


3 


4.6 X 10- 


3 


Z3w 


1.0 X 


10- 


2 


7 


3 X 


10- 


3 


4.6 X 


10- 


3 


2.2 X 


10- 


2 


4.4 


X 


10- 


3 


2 


X 


10- 


2 


3 


6 


X 


10- 


-3 


0.091 


5 


1 X 


10- 


3 


1.2 X 


10- 


2 


1.7 X 10- 


2 


Z5w 


2.5 X 


10- 


3 


1 


1 X 


10- 


3 


9.6 X 


10- 


4 


4.6 X 


10- 


3 


3.6 


X 


10- 


4 


3 


5 X 


10- 


3 


1 


1 


X 


10- 


-3 


0.054 


9 


5 X 


10- 


4 


1.5 X 


10- 


3 


2.5 X 10- 


3 


Z6w 


9.3 X 


10- 


4 


2 


X 


10- 


4 


5.8 X 


10- 


4 


1.7 X 


10- 


3 


5.3 


X 


10- 


5 


7 


7 X 


10- 


4 


7 


2 


X 


10- 


-4 


0.037 


2 


9 X 


10- 


4 


3.3 X 


10- 


4 


6.2 X 10- 


4 


Z6e 


7.0 X 


10- 


4 


1 


5 X 


10- 


4 


7.7 X 


10- 


5 


9.2 X 


10- 


4 


2.1 


X 


10- 


5 


4 


7 X 


10- 


4 


1 


5 


X 


10- 


-4 


0.036 


2 


1 X 


10- 


4 


1.2 X 


10- 


4 


3.3 X 10- 


4 


Z7e 


2.8 X 


10" 


4 


5 


3 X 


10- 


5 


3.4 X 


10- 


5 


3.7 X 


10- 


3 


4.3 


X 


10- 


6 


1 


3 X 


10- 


4 


1 


2 


X 


10- 


-4 


0.023 


6 


6 X 


10- 


5 


3.1 X 


10- 


5 


9.7 X 10- 


5 


Z5a 


5.4 X 


10" 


3 


2 


X 


10- 


3 


3.3 X 


10- 


3 


1.1 X 


10- 


2 


1.2 


X 


10- 


3 


5 


2 X 


10- 


3 


3 


3 


X 


10- 


-3 


0.070 


2 


4 X 


10- 


3 


4.1 X 


10- 


3 


6.4 X 10- 


3 


Z5b 


5.0 X 


10- 


3 


1 


8 X 


10- 


3 


3.3 X 


10- 


3 


1.0 X 


10- 


2 


1.4 


X 


10- 


3 


4 


9 X 


10- 


3 


3 


3 


X 


10- 


-3 


0.068 


2 


2 X 


10- 


3 


3.8 X 


10- 


3 


6.0 X 10- 


3 




Fig. 5. — The power spectrum density of the kinetic (sohd) and 
magnetic (dashed) energy densities, for two vertical net flux sim- 
ulations with Am = 1000 (run Zl, bold) and Am = 1 (run Z5, 
thin). Plotted are the shell integrated spectrum, represented by 
E'^ = 47rfc^|B(fc)p (where E denotes kinetic or magnetic energy 
density), normalized to background pressure Pq = Pnc?. The area 



poet 

enclosed by each curve corresponds to the total energy density from 
turbulent fluctuations. 



course discrete, but for notational convenience we write 
the formulas in continuous form. According to Parseval's 
theorem, we have 



- / IBix^d'x = 



\B{k)\ 



(27r)3 



(27r)3 



(22) 



Dividing by a factor of Stt we obtain the PSD for the 
magnetic energy density Mk = B'^/Stt. Similarly, one 
can obtain the PSD for the kinetic energy Ki^. 

In FigurelU show the PSDs computed from our runs 
Zl and Z5. These two simulations arc representative for 
the MRI turbulence in the ideal MHD and AD domi- 
nated regimes respectively, and arc run for two times 




Fig. 6. — Cumulative probability distribution of current density 
J for our simulations Zl {Am = 1000, black), Z3 {Am = 10, blue), 
Z5 [Am = 1, red). The current is normalized to •/^nPoc/AnH. 



longer than many other simulations (thus giving better 
statistics). We see that the shape of the PSD obtained 
from our simulations arc very similar. The PSD roughly 
follows a power-law form at small k, with the power law 
index approximately equals to —11/3, which is the in- 
dex for incompressible Kolmogorov turbulence spectrum. 
There appears to be a spectral break at kH 70, corre- 
sponding to a wavelength of about O.IH, and the PSD 
falls off rapidly toward smaller scales. The turbulent 
power in the Am — 1 case is about 20 times smaller than 
that in the ideal MHD case. Magnetic energy fluctu- 
ations dominate kinetic energy fluctuations in the ideal 
MHD case, while in the AD dominated regime, more tur- 
bulence power resides in the kinetic energy. Moreover, 
we have also checked the contour plot of vertically inte- 
grated PSD (not shown) and found that the turbulence 
becomes more anisotropic in the AD dominated regime: 
the turbulent power is more elongated in k^ than in ky. 

Second, we study the effect of AD on the distribution of 
current in the MRI turbulence. It has been shown that in 
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0.2 0.4 0.6 0.8 1 

X=|cos{J, B)| 



Fig. 7. — Cumulative probability distribution of |cosa|, where 
a is the angle between J and B. Shown are the results from our 
runs Zl {Am = 1000, black), Z2 {Am = 100, blue), Z3 {Am = 10, 
red). 

one and two dimensions, sharp current structure can be 
developed around magnetic nu lls in the presence of AD 
(|Brandenburg fc ZweibellHOOl . To examine whether the 
same effect is present in the MRI turbulence, we show in 
Figure |6] the cumulative probability distribution of the 
current density J = \ J\ in our simulation runs Zl, Z3 and 
Z5. If sharp current structure were to form, one would 
expect to see extended tails in the probability distribu- 
tion. However, we see that as Am decreases, the prob- 
ability distribution shifts leftward since turbulence be- 
comes weaker, but its shape remains largely unchanged. 
Wc also note that the current sharpening phenomenon is 
not observed in simulations by (|Brandenburg et al. 1993) 
either. It is likely the sharpening of current by AD is 
overwhelmed in 3D MHD turbulence. 

AD has also been shown to tend to reduce the cur- 
rent compone nt that is perpe ndicular to the magnetic 
field (Branden burg et al.lll995( ) and make the magnetic 
configuration more force-free. To examine this effect, we 
show the cumulative probability distribution of |cosq;| 
from our simulation runs Zl, Z2 and Z3 in Figure [7l 
where a is the angle between the current and the mag- 
netic field. The cumulative distribution functions from 
our runs Z4 and Z5 are almost identical to that from run 
Z3; where Am = 10. We confirm that AD makes the 
distribution more concentrated toward |cosa| = 1 (i.e., 
J 11 B)E Nevertheless, since the distribution of I cos a I 
is already peaked at 1 in the essentially ideal MHD run 
(Zl), the effect of AD does not modify the distribution 
of current orientation substantially. 

3.2. Net toroidal flux simulations 

^ We also re p ort a difference in our results from 
[Brandenburg et al.l l l 199511 . In their simulations, the distri- 
bution function peaks at |coso| = in ideal the MHD case, and 
AD concentrate the current toward | cos o| = 1. In our simulations, 
we find that the distribution function is already concentrated at 
I cos a I = 1 under ideal MHD, and AD simply makes it more 
concentrated. 



TABLE 3 
Net toroidal flux simulations. 



Run 


Am 


/3o 


Resolution 


Restart 


Turbule 


Yl 


1000 


100 


64 X 128 X 64 


No 


Yes 


Y2 


100 


100 


64 X 128 X 64 


No 


Yes 


Y3 


10 


100 


64 X 128 X 64 


Yl 


Yes 


Y4 


3.33 


100 


64 X 128 X 64 


Yl 


Yes 


Y5 


1 


100 


64 X 128 X 64 


Yl 


No 


Y6 


0.33 


100 


64 X 128 X 64 


Yl 


No 


Y7 


0.1 


400 


64 X 128 X 64 


Yl 


No 


Ylw 


1000 


100 


64 X 128 X 64 


No 


Yes 


Y2w 


100 


100 


64 X 128 X 64 


No 


Yes 


Y3w 


10 


400 


64 X 128 X 64 


Ylw 


Yes 


Y4w 


3.33 


400 


64 X 128 X 64 


Ylw 


Yes 


Y5w 


1 


400 


64 X 128 X 64 


Ylw 


No 


Y6w 


0.33 


400 


64 X 128 X 64 


Ylw 


No 


Y7w 


0.1 


400 


64 X 128 X 64 


Ylw 


No 


Ylh 


1000 


100 


128 X 256 X 128 


No 


Yes 


Y3h 


10 


100 


128 X 256 X 128 


Ylh 


Yes 


Y5h 


1 


100 


128 X 256 X 128 


Ylh 


No 


Y6h 


0.33 


100 


128 X 256 X 128 


Ylh 


No 



Box size is fixed at H X AH X H, u is fixed at 0.5. All simulations 
arc run for 24 orbits (ISOn^-*^). 

In the second group of simulations, the initial field 
configuration is chosen to be uniform along the az- 
imuthal direction y, with strength characterized by /3o = 
2Pn /Bn, where Bn is the in itial field strength. Follow- 
ing ISimon fc HawlevI (pOOl) . we fix the box size to be 
H X AH X H in the radial, azimuthal and vertical dimen- 
sions for all simulation runs in this groupQ- We choose 
the fiducial resolution to be 64 cells per H in the radial 
and vertical direction, and 32 cells per H in the azimuthal 
direction. All our net toroidal flux simulations are listed 
in Table [3l including one set of fiducial simulations with 
/3o = 100, one set of higher resolution simulations, and 
one set of weak field simulations with /3o = 400. Un- 
like net vertical flux, the net toroidal flux is not precisely 
co nserved in our shearing box simulations. As discussed 
in (jSimon fc HawlevI [20091) ■ ensuring strict conservation 
of toroidal flux numerically is more complex, and is also 
less important than conserving net vertical flux because 
the saturation level of the MRI turbulence is not very 
sensitive to the toroidal flux. Throughout all our simu- 
lations in this group, we find that the deviation of net 
toroidal flux from the initial value is generally less than 
2%. 

The linear stability of Keplerian disks in the pres- 
ence of pure toroidal field is more complex than that 
for the vertical field case. It requ ires consideration 
of no n-axisymmetric perturbations (jBalbus fc HawlevI 
Il992f ). and involves the time-dependent amplification of 
wave modes as the radial wave number swings from lead- 
ing to trailing. In ideal MHD, pure toroidal MRI fa- 
vors high kz wave numbers, and requires relatively large 
numerical resolution. In the case of Ohmic resistivity, 
swing amplification of modes is suppressed when the dif- 
fusion tim e of the mode is comparable to the orbital 
frequency ()Papaloizou fc Tergueml [19971 ). A linear sta- 

^ We have also performed our run series Yl to Y6 using a larger 
box 4H X AH X H and found that the turbulence properties are 
very similar. 
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Fig. 8. — The evolution of Maxwell stress in our fiducial set 
of toroidal net-flux runs (from top to bottom: Yl, Y2, Y7) 
normalized to CsH. For all runs after Y2, the simulations start 
from the end of the Yl run. 

bility with non-axisymmetric perturbations in AD domi- 
nated regime has yet to be performed. Nevertheless, one 
might expect that a similar argument holds for AD, with 
Am ^ 1 as the boundary for stability. 

3.2.1. A Fiducial Set of Runs 

We fix I3q = 100 and run 7 fiducial simulations with 
different Am values, named from Yl with Am — 1000 
to Y7 with Am = 0.1 (see Table [S]) similar to the case 
of net vertical flux runs. Initial growth of the MRI from 
pure toroidal field is more difficult and is only achieved 
when Am is greater than 10. We initialize the rest of the 
simulations from the turbulent state at the end of run 
Yl. 

Figure [5] illustrates the time evolution of the Maxwell 
stress from this fiducial set of runs. The time and volume 
averaged quantities from these runs are listed in Table 01 
Wc find that at a given value of Am, the saturation level 
of the MRI with net toroidal flux is much lower than 
the net vertical flux case. The turbulent energy density 
and total stress in run Yl (essencially ideal MHD) is a 
few times 10~^, about an order of magnitude less than 
run Zl. As Am drops below 10, the saturation level of 
the MRI turbulence falls off very rapidly. At Am ~ 3 
(run Y4), the total stress falls below 10"'^. At Am. ~ 1, 
although a total stress is maintained at a level of 10^'^, 
wc do not observe any signature of the MRI turbulence 
by examining the structure of the velocity field, which 
is essentially laminar (see further discussion in Section 
I3.2.4P . Since the simulation is initialized from a turbulent 
state, the low level of stress and kinetic energy are mostly 
due to the eigen-modes of the shearing box excited from 
the initial turbulence, and that are not damped due to 
the low dissipation in the Athena code. Unlike the case 
for net vertical fiux, there appears to exist a critical value 
of Am below which MRI turbulence with net toroidal flux 
is not self-sustained. This critical value of Am is about 
3. This fact will be further discussed shortly in Section 

HS also performed a number of net toroidal fiux two- 



fluid simulations with Am w 1 and Am w 100, and in 
both cases turbulence is self-sustained with total stress 
a on the order of 10"* to lO'^. A gain, these results 
are no longer valid in the strong-coupling regime and are 
not directly comparable to our results (see discussion in 
Section [HXT]) . 

We see from Table|3]that the density fluctuations in the 
net toroidal fiux simulations are generally smaller than 
those in the net vertical flux case. Therefore, following 
the discussion in Section 13.1.31 we expect the effect of 
ly has essentially no impact on the conclusions we have 
drawn above. 

3.2.2. The Effect of Resolution 

Relatively high resolution is needed for net toroidal 
flux MRI simulations in order to properly capture the 
amplificatio n of wave modes as the y swing from leading 
to trailing (jSimon fc Hawlevl 120091) . In order to justify 
our results in the previous subsection, we perform a few 
of the simulations with doubled resolution. These runs 
are labeled with an additional letter "h" (i.e., high reso- 
lution) in Tabled 

The time and volume averaged quantities from these 
high resolution runs are shown in Table HI We see 
that the kinetic and magnetic energy densities from in 
the high resolution simulations are generally larger than 
those in the low resolution runs, but only by a small 
factor. In particular, the difference between low and 
high resolutions with relatively large AD coefficient is 
only about 10% (e.g., comparing runs Y3 and Y3h), 
which strongly indicates numerical convergence. This is 
not surprising since small scale structures can be largely 
damped by AD thus higher resolution becomes unneces- 
sary. Run Y5h is also very similar to run Y5, where the 
initial turbulence is damped with remnant small velocity 
and magnetic fluctuations unlikely to be associated with 
the MRI turbulence. 

The inferences above are further justified by looking 
at the power spectrum of magnetic and kinetic energies. 
Following the same procedure described in Section [3 1.41 
we show in Figure |9] the shell integrated PSDs for runs 
Yl, Ylh and Y3, Y3h. The spectral shapes from the low 
resolution simulations appear to have a spectral peak 
at relatively large scales of about 0.5H, while at higher 
resolution, a power law spectrum at intermediate scales 
from approximately O.IH to 0.5-ff analogous to an inertia 
range appears to have developed. This observation may 
suggest high numerical resolution of at least 128 cells per 
H is needed for the toroidal field MRI simulations in or- 
der to resolve the inertia range in the turbulent spectrum, 
although smaller resolution of 64 cells per H appears to 
be sufficient for turbulence properties to converge. The 
shape of the PSDs at small k look very different from 
the PSDs in the net vertical fiux simulations, indicating 
different energy injection mechanism, while at large fc, 
the PSDs from both cases fall off in a similar manner, 
indicating similar dissipation mechanism. 

3.2.3. The Effect of toroidal Field Strength 

We have also performed the same set of net toroidal 
flux simulations with the toroidal magnetic field strength 
lowered by one half, /3q ~ 400, labeled with an additional 
letter "w" (i.e., weak field) in Tabic [3] General properties 
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TABLE 4 

Time and volume averaged quantities in net toroidal flux simulations. 



Run Ek,x Ek,y Ek,~ Ek Em,x l^M,y Em,z {Sp)/po QRc QMax " 



Yl 


1 


.2 


X 


10" 




1 


.1 


X 


10- 




5, 


.2 


X 


IQ- 




2, 


.8 


X 


10- 




9, 


.7 


X 


10- 


.3 


5, 


.0 


X 


10- 




4.4 


X 


IQ- 




0.10 


7, 


.4 


X 


10- 


.3 


2.6 


X 


10- 




3, 


.4 


X 


10- 


^27 


Y2 


1, 


.0 


X 


10" 


-2 


8, 


.2 


X 


10- 


-3 


4, 


.6 


X 


IQ- 


-3 


2, 


.3 


X 


10- 


-2 


8, 


.2 


X 


10- 


-3 


4, 


.0 


X 


10- 


-2 


4.0 


X 


10- 


-3 


0.083 


6, 


.2 


X 


10- 


-3 


2.0 


X 


10- 


-2 


2, 


.6 


X 


10- 


-2 


Y3 


4, 


.6 


X 


10- 


-3 


2, 


.4 


X 


10- 


'3 


1, 


.8 


X 


10- 


"3 


8, 


.8 


X 


10- 


-3 


2, 


.4 


X 


10- 


-3 


1, 


.9 


X 


10- 


-2 


1.4 


X 


10- 


-3 


0.062 


2, 


.5 


X 


10- 


-3 


5.7 


X 


10- 


-3 


8, 


.2 


X 


10- 


-3 


Y4 


6, 


.9 


X 


10" 


-4 


2 


.8 


X 


10- 


-4 


2, 


.1 


X 


IQ- 


-4 


1, 


.2 


X 


10- 


-3 


3, 


.4 


X 


10- 


-4 


1, 


.0 


X 


10- 


-2 


2.4 


X 


10- 


-4 


0.029 


3, 


.0 


X 


10- 


-4 


5.0 


X 


10- 


-4 


8, 


.1 


X 


10- 


-4 


Y5* 


3, 


.5 


X 


10" 


-5 


8 


.2 


X 


10- 


-5 


1, 


.8 


X 


IQ- 


-5 


1, 


.4 


X 


10- 


-4 


1, 


.4 


X 


10- 


-5 


9, 


.7 


X 


10- 


-3 


3.1 


X 


10- 


-5 


0.008 


6, 


.6 


X 


10- 


-6 


1.2 


X 


10- 


-5 


1, 


.9 


X 


10- 


-5 


Ylw 


7, 


.3 


X 


10" 


-3 


5 


.6 


X 


10- 


'3 


3, 


.1 


X 


IQ- 


-3 


1, 


.6 


X 


10- 


-2 


4, 


.8 


X 


10- 


-3 


2, 


.8 


X 


10- 


-2 


2.1 


X 


10- 


-3 


0.079 


4, 


.5 


X 


10- 


-3 


1.5 


X 


10- 


-2 


1, 


.9 


X 


10- 


-2 


Y2w 


4, 


.9 


X 


10- 


-3 


2 


.8 


X 


10- 


-3 


1, 


.8 


X 


IQ- 


-3 


9, 


.6 


X 


10- 


-3 


2, 


.5 


X 


10- 


-3 


1, 


.6 


X 


10- 


-2 


1.2 


X 


10- 


'3 


0.063 


2, 


.7 


X 


10- 


-3 


7.7 


X 


10- 


-3 


1, 


.0 


X 


10- 


-2 


Y3w 


1, 


.7 


X 


10- 


-3 


5, 


.8 


X 


10- 


-4 


4, 


.7 


X 


IQ- 


-4 


2, 


.7 


X 


10- 


-3 


4, 


.9 


X 


10- 


-4 


5, 


.4 


X 


10- 


-3 


2.6 


X 


10- 


-4 


0.043 


7, 


.2 


X 


10- 


-4 


1.5 


X 


10- 


-3 


2, 


.2 


X 


10- 


-3 


Y4w 


3, 


.3 


X 


10- 


-4 


8, 


.1 


X 


10- 


-5 


8, 


.1 


X 


IQ- 


-5 


4, 


.9 


X 


10- 


-4 


6, 


.8 


X 


10- 


-5 


3, 


.4 


X 


10- 


-3 


5.8 


X 


10- 


-5 


0.022 


1, 


.0 


X 


10- 


-4 


2.2 


X 


10- 


-4 


3, 


.2 


X 


10- 


-4 


Y5w* 


1, 


.2 


X 


10- 


-4 


6, 


.4 


X 


10- 


-5 


1, 


.4 


X 


IQ- 


-5 


2, 


.0 


X 


10- 


-4 


5, 


.6 


X 


10- 


-6 


2, 


.4 


X 


10- 


-3 


4.2 


X 


10- 


-6 


0.014 


2, 


.4 


X 


10- 


-5 


1.1 


X 


10- 


-5 


3, 


.5 


X 


10- 


-5 


Ylh 


1, 


.4 


X 


10- 


-2 


1, 


.7 


X 


10- 


-2 


6, 


.7 


X 


10- 


-3 


3, 


.7 


X 


10- 


-2 


1, 


.6 


X 


10- 


-2 


7, 


.3 


X 


10- 


-2 


8.1 


X 


10- 


'3 


0.12 


8, 


.4 


X 


10- 


-3 


3.8 


X 


10- 


-2 


4, 


.7 


X 


10- 


-2 


Y3h 


4, 


.5 


X 


10- 


-3 


2 


.9 


X 


10- 


-3 


1, 


.9 


X 


10- 


-3 


9, 


.3 


X 


10- 


-3 


3, 


.2 


X 


10- 


-3 


2, 


.3 


X 


10- 


'2 


1.9 


X 


10- 


'3 


0.056 


2, 


.5 


X 


10- 


-3 


6.8 


X 


10- 


-3 


9, 


.3 


X 


10- 


-3 


Y5h* 


7, 


.4 


X 


10- 


-5 


6 


.3 


X 


10- 


-5 


8, 


.7 


X 


10- 


-6 


1, 


.5 


X 


10- 


-4 


1, 


.0 


X 


10- 


-5 


1, 


.0 


X 


10- 


'2 


3.2 


X 


10- 


'5 


0.011 


1, 


.5 


X 


10- 


-5 


7.5 


X 


10- 


-6 


2, 


.2 


X 


10- 


'5 



*: These runs are not turbulent. 

from the saturated state of these runs are also shovifn 
in Table |4l We see that at the same value of Ato, the 
kinetic and magnetic energy density from the weak field 
simulations are smaller than those in our fiducial runs 
by a factor of 2-3. By inspection of the velocity field as 
well as the distribution of current density, we find that 
sustained turbulence is supported in run Y4w but not 
in run Y5w. Note that the time and volume averaged 
kinetic energy density from run Y4w is 4.9 x 10~*, which 
is slightly below our limit of 10~'^, but the total stress 
a = 3.2 X 10""* is reasonably large and is unlikely to be 
caused by inertia waves in the simulation box. Further 
discussion will be given in the next subsection 13.2.41 

Together with the results from Section I3.2.1[ we sum- 
marize the results from net toroidal flux simulations in 
Figure [TOl We see that for both values of /3o (100 and 
400), we do not observe any sustained MRI turbulence 
for Am < 3. Although we have not explored weaker net 
toroidal flux, based on the trend we see from the Fig- 
ure, even if MRI can be self-sustained at Am < 1 with 
weaker net toroidal flux, the resulting total stress a is 
unlikely to be above 10~^. This is in stark contrast with 
the net vertical flux simulations, and indicates that pure 
toroidal field geometry is more stable in the AD domi- 
nated regime. 

3.2.4. Cntenon for Sustained Turbulence 

A key objective of this paper is to study when MRI 
turbulence can be self-sustained in the presence of AD 
in the strong coupling limit. For the net vertical flux 
simulations, the criterion is relatively clear: turbulence 
can be sustained as long as the most MRI unstable mode 
fits into the simulation box. The situation is less clear for 
net toroidal flux simulations, and it remains to be seen 
whether the latter can be characterized as self-sustained 
turbulence. 

We take run Y5 as an illustrative example in this sub- 
section, but the analysis also applies to other marginal 
runs including Z6, Y5w, Y5h. In run Y5, after the initial 
turbulence has damped, the flow is largely laminar, ex- 
cept in a few very localized regions where some narrow 
but azimuthally elongated current structure is present 
and evolves very slowly with time. These features can 



be better demonstrated by computing the vertically in- 
tegrated Fourier power spectrum of magnetic and kinetic 
energies. Following the same procedure as described in 
Section 13.1.41 but integrating the full three-dimensional 
PSD over k^, we show the contour plot of the kinetic en- 
ergy PSD in the — ky plane from our runs Yl, Y4 and 
Y5 in Figure [m We see that in ideal MHD (Yl), the 
vertically integrated PSD has elliptic contours elongated 
and tilted toward the k^ axis. The contours are distorted 
and more elongated when AD is added (Y4). However, in 
run Y5, we see that the overall shape of the PSD contours 
are extremely elongated in the k^ direction. The origi- 
nal elliptic contours are almost destroyed, with irregular 
fragments distributed around the center. These irregular 
features in the vertically integrated PSD strongly indi- 
cate that the system is not in a turbulent state. We have 
also found that such irregular features are also present in 
runs Z6, Y5w and Y5h, but not present in runs Z5, Y4w 
and Y3h. 

In sum, we conclude that non-zero kinetic energy den- 
sity and total stress do not prove the existence of the 
MRI turbulence in shearing box simulations. Transition 
from self-sustained turbulence to non-turbulent state can 
be judged by looking at the vertically integrated PSD of 
kinetic and magnetic energies. 

3.3. Simulations with both vertical and toroidal fluxes 

Motivated by the linear stability analy sis of t h e MR I 
with AD by iKunz fc Balbug (|200l and iDeschI (|200l . 
we present simulations that include both vertical and 
toroidal fluxes. These authors found that in the pres- 
ence of both vertical and toroidal field, unstable modes 
exists for any values of Am, and the fastest growth rate 
is non- vanishing even as Am — >■ 0+. When Am < 1, 
the wave number of the most unstable mode has a sub- 
stantial non-zero radial component. In this subsection, 
we explore whether this behavior in the linear regime 
affects the non-linear evolution of the MRI with AD. 

We use the vertical plasma /?: /3zo = 8ttPo/BI to spec- 
ify the net vertical magnetic fiux. The net toroidal flux 
is specifled by its ratio to the net vertical fleld B^/Bz- 
We consider two sets of simulations, with B^/B^ = 4 
and B^/Bz = 1.25 and labeled by letters "M" and "N" 
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Fig. 9. — Similar to FigureO but for net toroidal flux simulations. 
Shown arc the shell integrated power spectrum density of the ki- 
netic (solid) and magnetic (dashed) energy densities. Upper panel: 
results from simulations with Am = 1000 (essentially ideal MHD). 
Our fiducial low-resolution results (run Yl) are shown in bold 
curves, while high-resolution results (run Ylh) are shown in thin 
lines. Lower panel: results from simulations with Am = 10, with 
low-resolution (run Y3) and high-resolution (Y3h) results shown 
in bold and thin lines. The uncertainty is about 10%. 



respectively. Parameters of these runs are given in Ta- 
ble [S] Similar to the previous subsections, we scan the 
parameter Am from 1000 do wn to 0.1. We use the dis- 
pcrsion relation (31) - (35) in lKunz fc BalbusI (|2004l) to 
find the wavenumbcr of the most unstable mode. We 
find that for Am > 3, the most unstable wavenumbcr is 
essentially purely vertical, while for Am < 1, the most 
unstable wavenumbcr is oblique with —kr. Corre- 

spondingly, we choose the box size to be i7 x AH x H for 
runs with Am > 3 and AH x AH x H for Am < 1. For 
simulations with Am < 1, we expect the fastest growth 
of the MRI to occur in the diagonal direction of the x — z 
plane. 

By default, we fix /3zo = 1600 in both sets of simu- 
lations. However, as Am falls below 1, wc find that the 
fastest growing mode will no longer fit into our simulation 
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Fig. 10. — The time and volume averaged total stress a from 
our net toroidal fiux simulations. Simulations with different net 
toroidal flux, characterized by the plasma /3o, are labeled by dif- 
ferent symbols and colors. The arrows for each /3o (as represented 
in red and blue) indicate the range of Am where MRI turbulence 
can be sustained. 
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Fig. 11. — Contour plot of the vertically integrated power spec- 
trum density of kinetic energy from our net toroidal flux simulation 
runs Yl (Am = 1000), Y4 {Am = 3) and Y5 {Am = 1). Neighbor- 
ing contours are separated by a factor of 10 in the power spectrum 
density. 
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TABLE 5 

Simulations with both net vertical and toroidal fluxes. 



Run 


Am 


Box size 


n 1 

Ao 




Orbi 


Ml 


1000 


H X 4H X H 


1600 


4 


48 


M2 


100 


H X AH X H 


1600 


4 


48 


M3 


10 


H X AH X H 


1600 


4 


36 


M4 


3.33 


H X AH X H 


1600 


4 


36 


M5 


1 


AH xAH X H 


1600 


4 


42 


M6 


0.33 


AH xAH X H 


8000 


4 


96 


M7 


0.1 


AH xAH X H 


30000 


4 


96 


Nl 


1000 


H X AH X H 


1600 


1.25 


48 


N2 


100 


H xAH X H 


1600 


1.25 


48 


N3 


10 


H xAH X H 


1600 


1.25 


48 


N4 


3.33 


H xAH X H 


1600 


1.25 


48 


N5 


1 


AH X AH X H 


1600 


1.25 


46 


N6 


0.33 


AH xAH X H 


1600 


1.25 


96 


N7 


0.1 


AH xAH X H 


8000 


1.25 


96 



The grid resolution is fixed at 64 cells per H. All simulations are 
initiated from seed perturbations, and generate sustained turbu- 
lence. 

^ : plasma fi from the vertical magnetic field. 
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box. Based on the results from Section 13.11 MRI turbu- 
lence would not be sustained in this case. Therefore, we 
increase to the value such that the most unstable 
mode just fit into the box. In the case of B^jB^, = 4, 
/3zo is increased to 8000 and 3 x 10^ for Am = 0.33 and 
Am = 0.1 respectively. For B^/B^ = 1.25, we increase 
/3zO to 8000 at Am = 0.1. Consequently, for all simula- 
tions in this subsection, MRI grows from the initial seed 
perturbations, and generates sustained turbulence after 
saturation. Below we discuss the properties of the MRI 
turbulence in the two sets of simulations. 

The first group of simulations (with B^/Bz — 4) arc 
run for at least 36 orbits. The initial growth of the MRI 
from runs M5 (Am = 1), M6 (Am = 0.33) and M7 
{Am = 0.1) are slower than the ideal MHD limit due 
to strong effect of AD: the fastest growth rates a^n pre- 
dicted from the linear dispersion relation for these simu- 
lation runs are 0.189f2-\ 0.1490"! and 0.136^^-^ respec- 
tively. This is to be compared with the case with pure 
vertical field with the same Am, where the correspond- 
ing am is 0.428rj-\ 0.2181}-^ and 0.07417-1 respectively 
The presence of both vertical and toroidal field gives a 
smaller growth rate at relatively large Am, but am de- 
creases much slower with decreasing Am than the pure 
vertical field case. At Am < 0.1, a field configuration 
with both net vertical and toroidal fluxes becomes more 
favorable than the pure net vertical field geometry by 
having substantially larger am- 

In the second group of runs, we choose B^/B^ = 1.25, 
which generates the fastest grow rate at Am < 1 com- 
pared with any other values. The fastest growth rates 
am for runs N5 {Am = 1), N6 {Am = 0.33) and NT 
{Am = 0.1) are 0.37117-\ 0.253n~^ and 0.2061}-^ re- 
spectively, which are nearly two times larger than those 
for B^/Bz = 4. 

Our simulations with Am < 1 (M5-M7, N5-N7) are 
particularly interesting because the fastest growing mode 
has a non-zero radial wave number \kr\ comparable to 
\kz\. During the linear growth stage, we observe axisym- 
mctric structures in the x — z plane similar to channel 



Fig. 12. — The distribution of current density in run M5 with 
Am = 1, B^/B, = 4 at the break down of the "channel flow" 
before saturating into turbulence (at time t = 570 ^'^). 

modes, but tilted toward the diagonal direction. These 
structures grow to a large amplitude, and finally break 
down into turbulence. As an example, we show in Fig- 
ure [T^] the distribution of current density right before the 
break down of the channel flow for run M5. In the turbu- 
lent state, one can still observe the emergence of struc- 
tures elongated in the diagonal direction in the x — z 
plane from time to time, which then fragment and in- 
ject kinetic energy into the system. For the simulations 
with Am < 0.33, these events lead to sporadic increase 
of kinetic energy and Reynolds stress on time scales of 
10 — 20 orbits. Due to such long time variability, we run 
these models for longer (to 96 orbits) and take the time 
average from about 56 orbits (35017-^) onward. 

The time and volume averaged properties of the MRI 
turbulence in all simulations with both net vertical and 
toroidal fluxes are listed in Table [51 In Figure [T^] we 
plot the total turbulent stress a as a function of Am 
from the two groups of simulations. We see that at 
Am > 1, where all simulations have the same net ver- 
tical flux Pzo — 1600, runs with B^/Bz = 4 generate 
slightly stronger MRI turbulence than simulations with 
Btj,/Bz = 1.25. This trend can be extrapolated down to 
zero net toroidal flux, as one can compare with results in 
runs Z3w and Z5w in Table O The dependence of tur- 
bulent strength on the toroidal flux is relatively weak, 
and when the net vertical flux in doubled, as in runs Zl 
to Z5, the strength of the turbulence becomes stronger 
than our corresponding runs Ml to M5. 

For simulations with Am < 1, we find that the turbu- 
lent stress a exceeds the maximum possible a attainable 
by the pure net vertical flux simulations. At Am = 0.1, 
the maximum value of a is 6.1 x 10"^, as compared to 
about 9.7x 10^^ from the pure net vertical flux case. This 
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TABLE 6 

Time and volume averaged quantities in simulations with both net vertical and net toroidal fluxes. 



Run Ek,x Ek,y Ek,z Ek Em,x Em.v Em,z {Sp)/po QRc QMax " 



Ml 


5 


6 


X 


10- 




5 


8 


X 


10- 


-Pi — 

z 


1 


9 


X 


10- 


-Pi — 

z 






0.13 




6 


7 


X 


10- 


-Pi — 

z 




0.30 




2 


8 


X 


10- 


-Pi — 

z 


0.26 


3 


4 


X 


10- 


-Pi — 

z 




0.17 






0.21 




M2 


4 


2 


X 


10" 


2 


3 


4 


X 


IQ- 


2 


1 


2 


X 


IQ- 


2 


8 


8 


X 


10- 


2 


4 


3 


X 


10- 


2 




0.17 




2 





X 


10- 


2 


0.16 


2 


4 


X 


10- 


2 




0.10 






0.13 




M3 


1 


6 


X 


10- 


2 


9 


2 


X 


10- 


3 


5 


7 


X 


10- 


3 


3 


1 


X 


10- 


2 


1 


1 


X 


10- 


2 


5 


3 X 10- 


2 


6 


3 


X 


10- 


3 


0.087 


9 


1 


X 


10- 


3 


2 


8 X 10- 


2 


3 


7 X 10- 


2 


M4 


9 


5 


X 


10" 


3 


4 


7 


X 


IQ- 


3 


3 


7 


X 


10- 


3 


1 


8 


X 


10- 


2 


3 


9 


X 


10- 


3 


2 


2 X 10- 


2 


3 


3 


X 


10- 


3 


0.075 


4 


7 


X 


10- 


3 


8 


7 X 10- 


3 


1 


3 X 10- 


2 


M5 


3 


9 


X 


10" 


3 


2 


8 


X 


10- 


3 


1 


4 


X 


10- 


3 


8 


1 


X 


10- 


3 


1 


3 


X 


10- 


3 


1 


3 X 10- 


2 


1 


7 


X 


10- 


3 


0.066 


1 


7 


X 


10- 


3 


2 


8 X 10- 


3 


4 


5 X 10- 


3 


M6 


1 


7 


X 


10" 


3 


2 


1 


X 


10- 


3 


2 


2 


X 


10- 


4 


4 





X 


10- 


3 


1 


3 


X 


10- 


4 


2 


6 X 10- 


3 


3 


1 


X 


10- 


4 


0.058 


5 


6 


X 


10- 


4 


4 


X 10- 


4 


9 


6 X 10- 


4 


M7 


1 


4 


X 


10- 


3 


4 





X 


10- 


4 


6 


5 


X 


10- 


5 


1 


9 


X 


10- 


3 


3 





X 


10- 


5 


7 


6 X 10- 


4 


7 


3 


X 


10- 


5 


0.050 


4 


8 


X 


10- 


4 


1 


3 X 10- 


4 


6 


1 X 10- 


4 


Nl 


3 


8 


X 


10- 


2 


3 


9 


X 


10- 


2 


1 


4 


X 


10- 


2 


9 


1 


X 


10- 


2 


4 


2 


X 


10- 


2 




0.20 




1 


8 


X 


10- 


2 


0.20 


2 


2 


X 


10- 


2 




0.12 






0.14 




N2 


2 


8 


X 


10- 


2 


2 


2 


X 


10- 


2 


9 


4 


X 


10- 


3 


5 


9 


X 


10- 


2 


2 


7 


X 


10- 


2 




0.11 




1 


3 


X 


10- 


2 


0.13 


1 


6 


X 


10- 


2 


6 


8 X 10- 


2 


8 


3 X 10- 


2 


N3 


1 


4 


X 


10- 


2 


6 


7 


X 


10- 


3 


4 


5 


X 


10- 


3 


2 


5 


X 


10- 


2 


7 


4 


X 


10- 


3 


3 


6 X 10- 


2 


4 


4 


X 


10- 


3 


0.080 


6 


9 


X 


10- 


3 


2 


X 10- 


2 


2 


7 X 10- 


2 


N4 


6 


3 


X 


10- 


3 


3 


4 


X 


10- 


3 


2 


9 


X 


10- 


3 


1 


3 


X 


10- 


2 


1 


7 


X 


10- 


3 


9 


2 X 10- 


3 


2 


1 


X 


10- 


3 


0.068 


2 


6 


X 


10- 


3 


4 


6 X 10- 


3 


7 


3 X 10- 


3 


N5 


3 


2 


X 


10- 


3 


1 


8 


X 


10- 


3 


1 


3 


X 


10- 


3 


6 


3 


X 


10- 


3 


5 


2 


X 


10- 


4 


4 


7 X 10- 


3 


1 


3 


X 


10- 


3 


0.065 


1 


2 


X 


10- 


3 


1 


8 X 10- 


3 


3 


X 10- 


3 


N6 


1 


2 


X 


10- 


3 


1 


1 


X 


10- 


3 


9 





X 


10- 


4 


3 


2 


X 


10- 


3 


1 


2 


X 


10- 


4 


1 


8 X 10- 


3 


8 


7 


X 


10- 


4 


0.057 


3 


8 


X 


10- 


4 


4 


8 X 10- 


4 


8 


6 X 10- 


4 


N7 


9 


9 


X 


10- 


4 


2 


4 


X 


10- 


4 


1 


1 


X 


10- 


4 


1 


3 


X 


10- 


3 


2 


9 


X 


10- 


5 


5 


7 X 10- 


4 


1 


8 


X 


10- 


4 


0.042 


3 





X 


10- 


4 


1 


6 X 10- 


4 


4 


7 X 10- 


4 
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Fig. 13. — The time and volume averaged total stress a from our 
simulations with both net vertical and net toroidal flux. Shown are 
results from two groups of simulations with different ratio B^/B~, 
are labeled by different symbols and colors. Dashed line repre- 
sents the maximum value of a obtained from pure net vertical flux 
simulations (Figure |4]l. 

is consistent with the hnear dispersion properties dis- 
cussed before: the presence of both vertical and toroidal 
field raises the maximum growth rate in the Am < 1 
regime. While the values of a given by the B^/B^ ~ 4 
group arc still larger than those in the B^/Bz = 1.25 
group, the latter group produces larger Maxwell stress. 
We note that for Am < 1, we have chosen the largest 
possible net flux such that the vertical extent of the sim- 
ulation box can fit only one most unstable mode. Ac- 
cording to the discussion in Section 13.1.21 the strength 
of the MRI turbulence from these simulations represents 
the highestpossible level at the given value of Am and 
the given magnetic field geometry for our box size, and 
may also approach the highest level in real disks. 

4. MRI WITH AMBIPOLAR DIFFUSION; DIAGNOSTICS 

In this section, we combine the results from all our 
simulations and further discuss the criteria for whether 
the MRI can be self-sustained with AD in a more general 



context. 

The MRI acts as a dynamo which amplifies initially 
weak fields. Although the saturati on mechanism of the 
MRI i s not well unders tood (but see lPessah fc GoodmanI 
(|2009[) : iPessahl (|2010[ )). the magnetic field energy at 
the saturated state scales with the net magnetic flux 
and i s generally below e q uipartition with t hermal en- 
ergy (jHawlev et all (|1995| ): iSano et al.l (|1998| ) as well as 
this paper). Given the field geometry, there should ex- 
ist a one-to-one correspondence between the initial field 
strength characterized by /3o and the final field strength 
characterized by (the space and time averaged gas to 
magnetic pressure at the saturated state), with (/3) < /3o 
for weak field due to the MRI dynamo, and gradually 
transiting to (/?) « Pq where the background field is too 
strong field to be destabilized. 

The quantity (/?) is also very useful for studying non- 
ideal MHD effects because the value it controls the 
relative importance of various no n-ideal MHD e ffects 
Ohm ic, HaU and AD, e.g., see IWardld ([2007[ ): [Bail 
2011| )). Henceforth, we shall consider (/3) as a main 
diagnostic quantity on the MRI turbulence. 

In Figure [HI we show the scattered plot of a and (/3) 
from all our simulations with sustained turbulence with 
different field geometries. We see that regardless of the 
initial field geometry and the value of Am, there is a re- 
markably tight correlation between the two quantities at 
the saturated state of the MRI turbulence. The correla- 
tion can be represented by 



This result is consistent with findings by (jHawlev et alj 
Il995| ) in ideal MHD simulations, and extends it to the 
non-ideal MHD regime. Analytical study of the sat- 
uration of the MRI with Ohmic resisti vity by para - 
sitic modes also predicts similar relations (jPessa 
More explicitly, this relation translates to 

W « iB^B^il + R) , (24) 

where R is the ratio of Reynolds to Maxwell stress (typi- 
cally w 1/3 in the ideal MHD case) . This relation implies 
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Fig. 14. — Scattered plot of the total turbulent stress a and the 
plasma /3 at the saturated state of the MRI turbulence from all our 
simulations. Simulations with different field geometries are marked 
by different symbols and colors as indicated in the legend, where 
B^/Bz = and B^/Bz = oo correspond to pure net vertical and 
pure net toroidal flux simulations respectively. Dashed line shows 
the fitting curve = \/2a. 
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Fig. 15. — Scattered plot of AD coefficient Am and the plasma 
/3 at the saturated state of the MRI turbulence from all our sim- 
ulations. Simulations with different field geometries are marked 
by different symbols and colors as indic ated in the legend. The 
dashed curve shows the fitting formula H25p as a lower bound of 
(/9> > /3min (Am). 



that the Maxwell stress is approximately a fixed fraction 
of the total magnetic energy, a reasonable result if the 
magnetic field is dominated by turbulent fluctuations. 
Only two points appear to deviate from this correlation, 
which correspond to net toroidal flux simulations with 
Am = 3. We see from Table |4] that in these two simu- 
lations, the magnetic energy is dominated by the back- 
ground toroidal field (i.e., in the transition where MRI 
is marginally sustained), therefore producing smaller 
than predicted. 

Next, we consider the relation between (/?) and Am. 
In Figure [15] we show scattered plot of Am and (/3) . We 
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Fig. 16. — Diagnostics of the MRI in the AD regime. At a give n 
Am, MRI is permitted when (/9) > /3niin(Am) (see equation I I25II ). 
and in the MRI permitted region, the stress a can be inferred given 
the field strength at the saturated state characterized by 

see that (/?) does not strongly correlate with Am, but 
also depends on the field geometry and the initial field 
strength. However, combining the simulation from all 
field geometries allow us to identify the lower bound of 
at a given Am, denoted by /3min, below which the 
field is too strong for to be destabilized based on our 
discussions before. For Am < 1, we have performed sim- 
ulations with the smallest possible /3o such that the most 
unstable mode marginally fit into the disk height H , and 
the value of /3,nin identified in this regime is robust. For 
Am > 1, the exploration on /3o may not be as complete 
especially in the simulations with both net vertical and 
toroidal fluxes and the actual /3min may be somewhat 
smaller than obtained here. Nevertheless, this regime is 
closer to ideal MHD and is less concerning. By com- 
bining all the available simulations, we obtain a fitting 
formula for /3niin given by 



/3niin(^TO) 



50 



Ai 



,1.2 



Ai 



,0.3 



1 



2-| 1/2 



(25) 



and is indicated in Figure [15] It asymptotes to 1 at 
Am — > cxj as one expects, while approaches 50/ Am^'^ for 
Am < 1. 

The constraint on /3nnn at a given Am allows us to 
identify the regions in the Am-{j3) plane at which MRI 
can or can not operate. In the mean time the correla- 
tion between a and (/3) provide the corresponding stress 
when MRI is permitted. Combining them together, the 
main results from the whole paper are best summarized 
in Figure 1161 MRI permitted regions are in the upper 
right with the boundary given by equation (j25p . It pro- 
vides useful diagnostics on the properties of the MRI in 
the AD regime in a concise fashion. 

First, at a given Am, the ultimate strength of the MRI 
turbulence (e.g., a and (/3)) depends on the field geome- 
try (including the net flux), but there exists a maximum 
a (or minimum (/?) ) at the most favorable field geometry 
(usually contains both net vertical and toroidal fluxes). 
One way to think about it is to starts with a weak reg- 



18 



Bai & Stone 



ular field as we perform our simulations. As the system 
evolves and as the MRI amplifies the field, the corre- 
sponding position of the system in the diagram moves 
downward and until it stops at some (/?) > /3min- 

Second, MRI can be self-sustained for any value of Am 
even for Am <C 1. Although we have explored the Am 
parameter down to Am = 0.1, we believe that it can be 
extended to further smaller Am because of the follow- 
ing reasons . Line ar analysis by iKunz &: Balbuj (|2004f ) 
and iDeschI (|2004[ ) shows in the presence of both verti- 
cal and toroidal field, MRI can grow at appreciable rate 
(approximately O.lSri""'^ when B^/Bz = 4) even in the 
limit of Am — )■ 0"*" provided that the field is sufficiently 
weak. This means that MRI turbulence can always be 
self-sustained. Meanwhile, we find that the linear dis- 
persion relation has already approached the small Am 
asymptote for Am < 0.3. Therefore, we expect the trend 
in Figure [15] on /3niin to hold to further smaller Am val- 
ues. 

Third, the boundary between the MRI permitted 
and prohibited regions is only suggestive but it does 
not necessarily imply sharp transitions. Our simula- 
tions are restricted by the limited box height (H) since 
they are unstratified. In reality, as one increases the 
field strength, the transition from sustained MRI tur- 
bulence to its suppression involves the effect of verti- 
cal stratification of gas density in the disks, and may 
be a smooth process. Before justified by stratified sim- 
ulations, which is left for our future work, this re- 
sult should be taken with some caution. In particu- 
lar, when vertical stratification is i nclude, linear analy- 
sis by lGammie fc BalbusI (|1994[ ) and lSalmeron &: Wardld 
(|2005D for ideal and non-ideal MHD have suggested the 
existence of global modes in the disk even in low /3o and 
small Elsasser number. On the other hand, in the case 
of Ohmic resistivity, the criterion that Ohmic Elsasser 
number equals one being the boundary between MRI 
permitted and su ppressed regions i dentified in unstrat- 
ified si mulations ( Sano et al.l 119981: [Fleming et all 120001: 
ISano fc Stouc 2002b|r do agree with resul t s from strat- 
ified simulations iFleming fc Stond [2001 [Turner et al.l 
[20?y7t lllgner fc Nelsonll20(m 

5. SUMMARY AND DISCUSSION 

Weakly ionized plasma is subject to a number of non- 
ideal MHD effects due to the coUisional coupling between 
the ionized species and the neutrals. Among them, am- 
bipolar diffusion (AD) results from the relative motion 
between the ions and the neutrals. It becomes most 
important when the gyro-frequencies of both the elec- 
trons and the ions in the magnetic field arc larger than 
their collision frequencies with the neutrals, so all ionized 
species are effectively coupled to the magnetic field. Con- 
sequently, AD usually dominates other non-ideal MHD 
effects (Ohmic resistivity and Hall effect) in regions with 
low density and high magnetic field. When the ion in- 
ertia is negligible and when the electron recombination 
time is much less than the dynamical time, AD is in the 
"strong coupling" limit, and can be studied by single- 
fiuid models. In this limit, the effect of AD is fully char- 
acterized by the parameter Am, the number of times a 
neutral molecule / atom collide with the ions in a dynam- 
ical time, and is the equivalence of the Elsasser number 



for Ohmic resistivity. The effect of AD becocmes dynam- 
ically important when Am approaches order unity. 

In weakly ionized disks such as the protoplanetary 
disks (PPDs), AD dominates other non-ideal MHD ef- 
fects in the d isk upp er layer as well as the outer regions 
of the PPDs ([WardlG,2007, : Bai.2011i ). The magnetorota- 
tional instability (MRI), which is considered as the major 
mechanism for providing angular momentum transport 
via the MHD turbulence, is strongly affected by AD. In 
the linear regime, the growth of the MRI is reduced or 
suppressed in the presence of AD, depen ding on the value 
of Am and the magnetic field geometry (|Blaes fc Balbu j 
[1991 IKunz fc Balbusl[200llDeschl[200l . Two-fluid sim- 
ula tions of the non-l i near e volution of the MRI with AD 
by [Hawlev fc Stond (|1998D showed that significant tur- 
bulence and angular momentum transport occurs when 
Am > 100. However, they ignored the processes of ion- 
ization and recombination, and the results are not di- 
rectly appli cable to th e PPDs, where the strong coupling 
limit holds ([Bai[[20lH ). 

We have implemented AD in the strong coupling limit 
in the Athena MHD code using a first-order accurate 
operator-split method. Its performance is verified by nu- 
merical tests on standing C-type shocks and the damping 
of MHD waves. We then perform local shearing box nu- 
merical simulations to study the effect of AD on the non- 
linear evolution of the MRI in the strong coupling limit. 
Our simulations are vertically unstratified with vertical 
box size equals the disk scale height H. The main pur- 
pose of this paper is to study how the properties of the 
MRI turbulence is affected by AD, especially the condi- 
tion under which MRI turbulence can be self-sustained. 
We perform three groups of simulations with different 
magnetic field configurations, and the main results are 
summarized below. 

1. Net vertical flux simulations. Unstable linear MRI 
modes always exist for any value of Am, with longer 
wavelength and smaller growth rate as Am gets 
smaller. MRI turbulence can be self-sustained as 
long as the wavelength of the most unstable mode 
Xm is within H (so the vertical field has to be pro- 
gressively weaker for Am — >■ 0+). At fixed Am, 
the total turbulent stress a increases monotonically 
with net vertical flux, until reaches the maximum 
when A,„ ~ H. The maximum value of a rapidly 
decreases with Am when Am < 10, from about 0.4 
at Am — )- oo down to about 0.007 at Am = 1. It 
falls below 10~^ at Am « 0.3 and is around of 10~^ 
at Am = 0.1. 

2. Net toroidal flux simulations. This field configu- 
ration is more stable. At fixed Am > 3 and net 
flux, the turbulent stress a from net toroidal flux 
simulations is smaller than that from net vertical 
flux simulations by about an order of magnitude. 
We do not flnd any evidence that MRI turbulence 
can be self-sustained at the level of a > 10~^ when 
Am < 1 for any net toroidal flux. 

3. Simulations with both net vertical and net toroidal 
fluxes. When Am > 1 and flxed net vertical flux, 
the strength of the MRI turbulence is similar to 
the pure net vertical flux case, but slowly increases 
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with the net toroidal flux. When Am < 1, the 
most unstable mode has non-zero radial wavenum- 
ber comparable to the vertical wavenumber, and 
the fastest growth rate asymptotes to some appre- 
ciable value even as Am — ^ 0"*". The maximum 
value of the turbulent stress a largely exceeds the 
pure net vertical flux case when Am < 1, with 
a « 6 X 10"'' at Am = 0.1. 

In addition, we find that similar to the effect of Ohmic 
dissipation, the ratio of the fluctuating part of the mag- 
netic energy density to the kinetic energy density de- 
creases as AD is stronger. Similarly, the ratio of Maxwell 
stress to Reynolds stress also drops at smaller Am. The 
power spectra density of the MRI turbulence in the AD 
dominated regime does not show any new features other 
than a rescaling from that in the ideal MHD case. We 
do not find any evidence that AD leads to the forma- 
tion of sharp c urrent structures in the MRI tu rbulence 
as proposed by [Brandenburg fc Zweibel (jl994[ ). but we 
confirm that AD tends to reduce the component of cur- 
rent p erpendi cular to th e direction of magnetic field 
(jBrand cnburg ct al.lll995[) . although to a lesser extent. 

Combining the results from these three groups of sim- 
ulations, we find a strong correlation between the tur- 
bulent stress a and the gas to magnetic pressure ratio (3 
at the saturated state, given by {/3) w l/2a. The sus- 
tainability and saturation level of the MRI turbulence 
depends on the value of Am, the magnetic field geome- 
try, and the magnetic field strength. It is best summa- 
rized in Figure [TBI In short, at a given Am, there exists a 
maximum value of turbulent stress a achievable from the 
most favorable geometries (generally with both net verti- 
cal and toroidal fluxes). Correspondingly, at a given Am, 
there exists a maximum field strength above which MRI 
is suppressed, and the maximum field strength rapidly 
decreases with decreasing Am. For future reference, we 
quote the turbulent stress a = 7 x 10""^ and 6 x 10"'' as 
the maximum value we have found in our simulations at 
Am = 1 and 0.1 respectively. 

In principle, in the presence of both net vertical and 
net toroidal fluxes, MRI turbulence can be self-sustained 
at any values of Am, provided that the magnetic field 
is sufficiently small. Nevertheless, the resulting turbu- 
lent stress a would be much smaller than 6 x 10~^ when 
Am < 0.1, and would be inefficient in transporting an- 
gular momentum in most astrophysical disks. Moreover, 
we recall that AD dominates other non-ideal MHD cf- 
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fects (Ohmic resistivity and Hall effect) at relatively low 
density and relatively large magnetic field. When the 
magnetic field is too weak as required for the wavelength 
of the most unstable MRI mode to be within H, AD may 
no longer be the dominant non-ideal MHD effect, and our 
results are not directly applicable. Therefore, while gen- 
eralization of our results to Am < 0.1 are possible, it 
may not be physically relevant. 

We have shown that the effect of AD on the non-linear 
evolution of the MRI depends on the field geometry, 
which is arbitrarily assigned in our local shearing box 
simulations. In reality, the field geometry depends on 
the global evolution of the disk. Let us imagine one sce- 
nario and apply our local simulation results to a global 
picture. If the disk is initially thre aded by a very weak 
vertical field (e.g., lArmitagg 119981 ). the initial growth 
of the MRI may be governed by the Ohmic resistivity 
and/or the Hall effect. AD takes over when substantial 
amplification of the field occurs and dominate the non- 
linear evolution of the MRI. The saturation of the MRI 
generates strong vertical and azimuthal fields, and redis- 
tributes the vertical and toroidal fluxes across the entire 
disk. No longer limited by the enforced net fluxes as in 
our local simulations, the saturation of the MRI could 
probably arrange the magnetic field to achieve the most 
favorable local field configurations such that the field is 
maximumly amplified and (/3) « /3min is achieved. There- 
fore, given the value of Am in the disk, and provided that 
AD dominates other non-ideal MHD effects, we expect 
the field strength in the disk to be largely determined by 
the value of PminiAm), via Equation (j25p . from which 
one can also obtain a rough estimate of a « l/2/3min- 

Our results are mostly relevant to the structure and 
evolution of the PPDs. Details about the application 
require considerations of the ionization and recombina- 
tion processes in the disks with an appropriate chemistry 
model, which is be yond the scope of this paper. In our 
companion paper (iBaill201lD . we will take into account 
all non-ideal MHD effects together and study their im- 
plications in the PPDs. 
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